Data Input and Manipulation

# update.packages(ask = FALSE, checkBuilt = TRUE)  # update R packages
# tinytex::install_tinytex()
# tinytex::tlmgr_update()  # update LaTeX packages
#tinytex::reinstall_tinytex(force = TRUE)
# setwd("~/Dropbox/PhD/PAPERS/3_n1_tim/sub.EE/timema-eco-evo-2020/")
library(plotrix)
## Warning: package 'plotrix' was built under R version 3.6.2
# all arthropod data 
ad <- read.csv( "data/all_art_n1_data2013.csv", colClasses=c(plant_id = "factor"))
ad <- droplevels(ad[ad$plant_id != "021", ])

# plant-level data
n1 <- read.csv("data/n1_plant_data_16Sept14.csv", colClasses=c(plant_id="character"))
n1$plant_id <- as.factor(ifelse(nchar(n1$plant_id) == 1, paste("00", n1$plant_id, sep = ""), 
                      ifelse(nchar(n1$plant_id) == 2, paste("0", n1$plant_id, sep=""), n1$plant_id)))

# carbon-nitrogen data
cnd.raw <- read.csv("data/cn_data.csv")
cnd.raw$id <- as.factor(ifelse(nchar(cnd.raw$id) == 1, paste("00", cnd.raw$id, sep = ""), 
                      ifelse(nchar(cnd.raw$id) == 2, paste("0", cnd.raw$id, sep=""), cnd.raw$id)))
nfmal <- function(mal=NULL, k=1) {
  freq <- 1 - mal
  # adapted moph
  min.RSS <- function(data, par) {
              with(data, sum(( (par[1] / (1 - ((par[2] - par[1])/par[2]) * exp((-1) * k * x))) - y)^2))
  }

  dat.ad <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
  pars.ad <- optim(par = c(0,1), min.RSS, data=dat.ad)$par
  
  freqs <- seq(0, 1, .01)
  ws.ad <-  (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freqs)))
  
  # plot(freqs, ws.ad)
  
  # maladapted morph
  
  ws.mal <- rev(ws.ad) - .35 
  
  raw <- (ws.ad * freqs) + (ws.mal * (1 - freqs))
  
  min.w <- min(raw)
  max.w <- max(raw - min.w)
  
  scaled <- (raw - min.w)/max.w
  cbind(freqs, scaled)
  plot(rev(freqs), 1 - scaled)
  #plot(freqs, ws.ad, ylim=c(-2, 2), type="l")
  #plot(freqs, ws.ad, type="l")
  #lines(freqs, ws.mal, lty=2)
  
  # for output
   ws.ad2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freq)))
   ws.mal2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * (1-freq)))) - .35
   raw2 <- (ws.ad2 * freq) + (ws.mal2 * (1 - freq))
   sc2 <- (raw2-min.w) / max.w
  
  #plot(freqs, mal)
  #pars.ad

  1 - sc2

}
## overwrite with new maladaptation measure

# manipulate data and add variables 
n1$cn2 <- cnd.raw$cn[match(n1$plant_id, cnd.raw$id)]
n1$pC <- cnd.raw$wC[match(n1$plant_id, cnd.raw$id)]
n1$pN <- cnd.raw$wN[match(n1$plant_id, cnd.raw$id)]

n1$host <- factor(n1$host)
n1$occupied <- ifelse(n1$no_tim == 0, 0, 1)
n1$pcon_mal <- ifelse(n1$host=="A", n1$con_tim_C/n1$con_tim, n1$con_tim_A/n1$con_tim)
n1$con_mal <- ifelse(n1$host=="A", n1$con_tim_C, n1$con_tim_A)
n1$con_ad <- ifelse(n1$host=="C", n1$con_tim_C, n1$con_tim_A)
n1$c.con_art5 <- n1$con_art5-mean(n1$con_art5)
n1$c.con_art <- n1$con_art-mean(n1$con_art)
n1$nf_mal_k0.01 <- nfmal(n1$mal, k=.01)

n1$nf_mal_k0.6 <- nfmal(n1$mal, k=0.6)

n1$nf_mal_k0.5 <- nfmal(n1$mal, k=0.5)

n1$nf_mal_k1.2 <- nfmal(n1$mal, k=1.2)

n1$nf_mal_k2.0 <- nfmal(n1$mal, k=2)

n1$nf_mal_k3.0 <- nfmal(n1$mal, k=3)

n1$nf_mal_k4.0 <- nfmal(n1$mal, k=4)

n1$c.mal <- n1$mal-mean(n1$mal, na.rm=TRUE)
n1$c.nf_mal_k0.01 <- n1$nf_mal_k0.01 - mean(n1$nf_mal_k0.01, na.rm=TRUE)
n1$c.nf_mal_k0.6 <- n1$nf_mal_k0.6 - mean(n1$nf_mal_k0.6, na.rm=TRUE)
n1$c.nf_mal_k0.5 <- n1$nf_mal_k0.5 - mean(n1$nf_mal_k0.5, na.rm=TRUE)
n1$c.nf_mal_k1.2 <- n1$nf_mal_k1.2 - mean(n1$nf_mal_k1.2, na.rm=TRUE)
n1$c.nf_mal_k2.0 <- n1$nf_mal_k2.0 - mean(n1$nf_mal_k2.0, na.rm=TRUE)
n1$c.nf_mal_k3.0 <- n1$nf_mal_k3.0 - mean(n1$nf_mal_k3.0, na.rm=TRUE)
n1$c.nf_mal_k4.0 <- n1$nf_mal_k4.0 - mean(n1$nf_mal_k4.0, na.rm=TRUE)
n1$no_mal <- ifelse(n1$host=="C", n1$no_str_tim, n1$no_grn_tim)
n1$no_adp <- ifelse(n1$host=="A", n1$no_str_tim, n1$no_grn_tim)
n1$con_grn <- n1$con_grn_A + n1$con_grn_C 
n1$con_str <- n1$con_str_A + n1$con_str_C
n1$rcon_grn <- n1$con_grn/n1$con_tim
n1$rcon_grn_A <- n1$con_grn_A/n1$con_grn
n1$rcon_grn_C <- n1$con_grn_C/n1$con_grn
n1$con_q2N_A <- n1$q2N_A*n1$con_tim
n1$con_q2N_C <- n1$q2N_C*n1$con_tim
n1$rcon_q2N_C <- n1$q2N_C/(n1$q2N_A+n1$q2N_C)
n1$rcon_vol_C <- n1$con_vol_C/n1$con_vol
n1$rcon_tim_C <- n1$con_tim_C/n1$con_tim
n1_occ <- n1[n1$no_tim_gs>0,]
n1$dens_art_5 <- n1$no_art_5 / n1$ln.vol_cub
n1$dens_tim_5 <- n1$no_tim_5 / n1$ln.vol_cub
n1$cn_rat <- n1$cn_ratio
n1$vol_cubm <- n1$vol_cub * 0.000016
n1$any_dens <- n1$no_any/n1$vol_cubm
n1$any_dens_4 <- n1$no_any_4/n1$vol_cubm
n1$mal2 <- n1$mal^2
n1$c.mal2 <- n1$c.mal^2
n1$c.no_art_5 <- n1$no_art_5 - mean(n1$no_art_5)
n1$c.no_tim <- n1$no_tim - mean(n1$no_tim)
n1$c.host <- ifelse(n1$host == "A", -.5, .5)
n1$con_q2N_mal <- ifelse(n1$host == "A", n1$con_q2N_C, n1$con_q2N_A)
n1$con_q2N_ad <- ifelse(n1$host == "C", n1$con_q2N_C, n1$con_q2N_A)
# subset data -- using arthropods over 5mm in body length only
dd_5 <- droplevels(subset(ad, !is.na(ad$weight) & ad$is_5 ))
dd_5$plant_id <- as.factor(dd_5$plant_id)


bins <- 20 # choose number of bins
dd_5$mass.bin<-cut(log(dd_5$weight), bins) # cut into pieces
work <- dd_5[, c("plant_id", "host", "weight", "mass.bin")] # extract columns


pids <- levels(work$plant_id)
ll <- length(pids)

# initialize structures to collect coefficients
slopes<-rep(NA, ll)
int1 <- rep(NA, ll)
int3 <- rep(NA, ll)
int5 <- rep(NA, ll)
int7 <- rep(NA, ll)
int8 <- rep(NA, ll)
int9 <- rep(NA, ll)
int10 <- rep(NA, ll)
int18 <- rep(NA, ll)
int20 <- rep(NA, ll)
pid <- levels(work$plant_id)
sd <- data.frame(pid, int1, int3, int5, int7, int8, int9, int10, int18, int20, slopes)

# loop over plants to create mass-abundance relationships
for(i in 1:ll) {
  
  plant <- subset(work, plant_id==pids[i])

  # get counts of artropods in each bin
  mass.ab <- data.frame(with(plant, table(mass.bin)), cat=levels(plant$mass.bin))
  
  # create rank of log mass bins, and recenter data for intercept analysis
  mass.ab$cat.code0 <- 1:bins
  mass.ab$cat.code1 <- mass.ab$cat.code0 - 1 
  mass.ab$cat.code3 <- mass.ab$cat.code0 - 3
  mass.ab$cat.code5 <- mass.ab$cat.code0 - 5
  mass.ab$cat.code7 <- mass.ab$cat.code0 - 7
  mass.ab$cat.code8 <- mass.ab$cat.code0 - 8
  mass.ab$cat.code9 <- mass.ab$cat.code0 - 9
  mass.ab$cat.code10 <- mass.ab$cat.code0 - 10
  mass.ab$cat.code18 <- mass.ab$cat.code0 - 18
  mass.ab$cat.code20 <- mass.ab$cat.code0 - 20
  
  # replace 0 with NAs
  mass.ab$Freq[mass.ab$Freq==0] <- NA
  plot.these <- mass.ab
  
  # run simple regression models
  mod0 <- lm(log(Freq) ~ cat.code0, data = plot.these)
  mod1 <- lm(log(Freq) ~ cat.code1, data = plot.these)
  mod3 <- lm(log(Freq) ~ cat.code3, data = plot.these)
  mod5 <- lm(log(Freq) ~ cat.code5, data = plot.these)
  mod7 <- lm(log(Freq) ~ cat.code7, data = plot.these)
  mod8 <- lm(log(Freq) ~ cat.code8, data = plot.these)
  mod9 <- lm(log(Freq) ~ cat.code9, data = plot.these)
  mod10 <- lm(log(Freq) ~ cat.code10, data = plot.these)
  mod18 <- lm(log(Freq) ~ cat.code18, data = plot.these)
  mod20 <- lm(log(Freq) ~ cat.code20, data = plot.these)
  
  # collect the coefficients
  sd$int1[i] <- coef(mod1)[1]
  sd$int3[i] <- coef(mod3)[1]
  sd$int5[i] <- coef(mod5)[1]
  sd$int7[i] <- coef(mod7)[1]
  sd$int8[i] <- coef(mod8)[1]
  sd$int9[i] <- coef(mod9)[1]
  sd$int10[i] <- coef(mod10)[1]
  sd$int18[i] <- coef(mod18)[1]
  sd$int20[i] <- coef(mod20)[1]
  sd$slopes[i]<-coef(mod0)[2]
  
}

# fix plant ID codes to be 3 character factor
sd$pid <- as.character(sd$pid)
sd$pid <- as.factor(ifelse(nchar(sd$pid) == 1, paste("00", sd$pid, sep = ""), 
                                ifelse(nchar(sd$pid) == 2, paste("0", sd$pid, sep=""), 
                                       sd$pid)))

# NA all intercepts for which no slope could be calculated (n = 1 or 2)
sd[is.na(sd$slopes), 2:11] <- NA

# add data to n1

## add intercepts to n1

n1$mas <- NA
n1$mai1 <- NA
n1$mai3 <- NA
n1$mai5 <- NA
n1$mai7 <- NA
n1$mai8 <- NA
n1$mai9 <- NA
n1$mai10 <- NA
n1$mai18 <- NA
n1$mai20 <- NA

# add data, excluding plant 021
n1[match(sd$pid[sd$pid != "021"], n1$plant_id), 
   c("mai1", "mai3", "mai5", "mai7", "mai8", "mai9", "mai10", "mai18", "mai20", "mas")] <- sd[sd$pid != "021", 2:11]
bins <- 20 # choose number of bins
dd_5$mass.bin<-cut(log(dd_5$weight), bins) # cut into pieces

work2 <- dd_5[, c("plant_id", "host", "weight", "mass.bin")]

mass.ab <- data.frame(table(dd_5$mass.bin))

# check it
sum(mass.ab$Freq)
## [1] 1102
plot(x=1:20, y=mass.ab$Freq)

mad <- droplevels(ad[!is.na(ad$weight) & ad$is_4, ])

mas <- aggregate(mad$weight, by=list(plant_id = mad$plant_id), FUN=function(x) {
  
  if(length(x) < 3) return(NA)
  
  lw <- log(x)
  h1 <- hist(x)
  cnt <- h1$counts
  cnt[cnt==0] <- NA
  mids <- h1$mids
  an1 <- lm(log(cnt) ~ log(mids))
 
  return(c(mas=an1$coefficients[2]))
  
})

# add to N1 data
n1$mas <- mas[match(n1$plant_id, mas$plant_id), "x"]

Regression Analyses

Drivers of Maladaptation

an1 <- lm(mal ~ host, data=n1)
summary(an1)
## 
## Call:
## lm(formula = mal ~ host, data = n1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53932 -0.27745 -0.03932  0.26068  0.72255 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.27745    0.06669   4.160  7.2e-05 ***
## hostC        0.26187    0.08040   3.257  0.00158 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3592 on 91 degrees of freedom
##   (53 observations deleted due to missingness)
## Multiple R-squared:  0.1044, Adjusted R-squared:  0.09457 
## F-statistic: 10.61 on 1 and 91 DF,  p-value: 0.001582

Effects of Maladaptation

contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)

# main effect models
an1.1 <- glm(no_tim ~ c.mal + host + ln.vol_cub + con_tim + cn2,
             family=quasipoisson, 
           data=n1, subset=no_tim > 4)
summary(an1.1)
## 
## Call:
## glm(formula = no_tim ~ c.mal + host + ln.vol_cub + con_tim + 
##     cn2, family = quasipoisson, data = n1, subset = no_tim > 
##     4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.49224  -0.83543  -0.09667   0.62544   2.36662  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.207594   1.213446  -0.171   0.8652  
## c.mal       -0.028035   0.298378  -0.094   0.9257  
## host1       -0.045736   0.358165  -0.128   0.8992  
## ln.vol_cub   0.167263   0.073955   2.262   0.0304 *
## con_tim      0.051004   0.028218   1.807   0.0798 .
## cn2          0.001773   0.035684   0.050   0.9607  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.90865)
## 
##     Null deviance: 35.583  on 38  degrees of freedom
## Residual deviance: 29.055  on 33  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
an1.2 <- glm(no_tim ~ c.mal + host + ln.vol_cub + con_tim,
             family=quasipoisson, 
           data=n1, subset=no_tim > 4)
summary(an1.2)
## 
## Call:
## glm(formula = no_tim ~ c.mal + host + ln.vol_cub + con_tim, family = quasipoisson, 
##     data = n1, subset = no_tim > 4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.51373  -0.80168  -0.06667   0.59580   2.36304  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.149383   0.952609  -0.157   0.8763  
## c.mal       -0.005711   0.286871  -0.020   0.9842  
## host1       -0.042412   0.154675  -0.274   0.7855  
## ln.vol_cub   0.164847   0.071695   2.299   0.0276 *
## con_tim      0.056262   0.024962   2.254   0.0306 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.8610003)
## 
##     Null deviance: 36.470  on 39  degrees of freedom
## Residual deviance: 29.299  on 35  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
an1.3 <- glm(no_tim ~ host + ln.vol_cub + con_tim,
             family=quasipoisson, 
           data=n1, subset=no_tim > 4)
summary(an1.3)
## 
## Call:
## glm(formula = no_tim ~ host + ln.vol_cub + con_tim, family = quasipoisson, 
##     data = n1, subset = no_tim > 4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.51720  -0.80017  -0.06795   0.59200   2.36087  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.14857    0.93844  -0.158   0.8751  
## host1       -0.04130    0.14220  -0.290   0.7731  
## ln.vol_cub   0.16473    0.07046   2.338   0.0251 *
## con_tim      0.05629    0.02458   2.290   0.0280 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.8370613)
## 
##     Null deviance: 36.47  on 39  degrees of freedom
## Residual deviance: 29.30  on 36  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# final model
an1.f <- glm(no_tim ~ ln.vol_cub + con_tim,
             family=quasipoisson, 
           data=n1, subset=no_tim > 4)
summary(an1.f)
## 
## Call:
## glm(formula = no_tim ~ ln.vol_cub + con_tim, family = quasipoisson, 
##     data = n1, subset = no_tim > 4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.51357  -0.79264  -0.04083   0.53831   2.37497  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.23088    0.88708  -0.260   0.7961  
## ln.vol_cub   0.17084    0.06667   2.563   0.0146 *
## con_tim      0.05505    0.02397   2.297   0.0274 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.8169915)
## 
##     Null deviance: 36.470  on 39  degrees of freedom
## Residual deviance: 29.371  on 37  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
##### With Neg-Freq Dep Mal Calc -- no difference

an1.1 <- glm(no_tim ~  con_tim + ln.vol_cub,
             family=quasipoisson, 
           data=n1, subset=no_tim > 4)
summary(an1.1)
## 
## Call:
## glm(formula = no_tim ~ con_tim + ln.vol_cub, family = quasipoisson, 
##     data = n1, subset = no_tim > 4)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.51357  -0.79264  -0.04083   0.53831   2.37497  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.23088    0.88708  -0.260   0.7961  
## con_tim      0.05505    0.02397   2.297   0.0274 *
## ln.vol_cub   0.17084    0.06667   2.563   0.0146 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.8169915)
## 
##     Null deviance: 36.470  on 39  degrees of freedom
## Residual deviance: 29.371  on 37  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)

#### for new maladaptation model -- from p = 0.03 to 0.056
an2.1 <- glm(no_art_5 ~ c.nf_mal_k0.5 + ln.vol_cub + c.con_art5*host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
## 
## Call:
## glm(formula = no_art_5 ~ c.nf_mal_k0.5 + ln.vol_cub + c.con_art5 * 
##     host, family = quasipoisson, data = n1, subset = no_tim > 
##     0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5517  -1.6891  -0.7854   0.7949   8.0273  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.39575    1.07406  -3.162  0.00216 ** 
## c.nf_mal_k0.5     0.40763    0.23405   1.742  0.08511 .  
## ln.vol_cub        0.46044    0.08406   5.478 4.13e-07 ***
## c.con_art5        0.06330    0.02624   2.413  0.01793 *  
## host1             0.41085    0.21784   1.886  0.06262 .  
## c.con_art5:host1 -0.05911    0.03704  -1.596  0.11416    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.568432)
## 
##     Null deviance: 652.94  on 92  degrees of freedom
## Residual deviance: 414.74  on 87  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
####


an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host, 
##     family = quasipoisson, data = n1, subset = no_tim > 0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.6751  -1.6447  -0.7855   0.8684   7.8261  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.07216    1.02646  -2.993   0.0036 ** 
## c.mal             0.49991    0.23689   2.110   0.0377 *  
## ln.vol_cub        0.43164    0.08042   5.367 6.55e-07 ***
## c.con_art5        0.06731    0.02611   2.578   0.0116 *  
## host1             0.47820    0.21732   2.200   0.0304 *  
## c.con_art5:host1 -0.06740    0.03701  -1.821   0.0720 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
## 
##     Null deviance: 652.94  on 92  degrees of freedom
## Residual deviance: 407.37  on 87  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an2.2 <- glm(no_art_5 ~ c.mal + ln.vol_cub + host * c.con_art5, 
             family=quasipoisson, data=n1, subset=no_tim>1)
summary(an2.2)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + host * c.con_art5, 
##     family = quasipoisson, data = n1, subset = no_tim > 1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8730  -1.9463  -0.5069   0.9554   7.5987  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      -1.83738    1.27432  -1.442  0.15454   
## c.mal             0.44057    0.33068   1.332  0.18780   
## ln.vol_cub        0.33915    0.09909   3.423  0.00112 **
## host1             0.19255    0.35264   0.546  0.58707   
## c.con_art5        0.05499    0.02912   1.888  0.06382 . 
## host1:c.con_art5 -0.03657    0.05358  -0.683  0.49755   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 6.056629)
## 
##     Null deviance: 431.96  on 65  degrees of freedom
## Residual deviance: 324.94  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an2.3 <- glm(no_art_5 ~ mal + host + ln.vol_cub + host*c.con_art5 + cn_rat, 
             family=quasipoisson, data=n1)
summary(an2.3)
## 
## Call:
## glm(formula = no_art_5 ~ mal + host + ln.vol_cub + host * c.con_art5 + 
##     cn_rat, family = quasipoisson, data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8894  -1.5379  -0.5226   0.9420   7.5316  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.69695    1.49892  -1.799  0.07580 .  
## mal               0.57739    0.25635   2.252  0.02708 *  
## host1             0.76222    0.40843   1.866  0.06572 .  
## ln.vol_cub        0.45973    0.08489   5.416 6.39e-07 ***
## c.con_art5        0.08046    0.02808   2.865  0.00534 ** 
## cn_rat           -0.04820    0.04701  -1.025  0.30829    
## host1:c.con_art5 -0.06674    0.04068  -1.641  0.10484    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.285246)
## 
##     Null deviance: 624.76  on 85  degrees of freedom
## Residual deviance: 361.47  on 79  degrees of freedom
##   (60 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an2.4 <- glm(no_art_5 ~ mal * host + ln.vol_cub, family=quasipoisson, data=n1)
summary(an2.4)
## 
## Call:
## glm(formula = no_art_5 ~ mal * host + ln.vol_cub, family = quasipoisson, 
##     data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.1205  -1.7340  -0.8982   0.8801   7.5285  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.13525    1.00302  -3.126   0.0024 ** 
## mal          0.51650    0.31353   1.647   0.1031    
## host1        0.64591    0.29033   2.225   0.0287 *  
## ln.vol_cub   0.41083    0.07754   5.298 8.54e-07 ***
## mal:host1   -0.29777    0.49148  -0.606   0.5462    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.511588)
## 
##     Null deviance: 652.94  on 92  degrees of freedom
## Residual deviance: 437.29  on 88  degrees of freedom
##   (53 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# main effect models
an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + cn2 + c.con_art5*host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + cn2 + c.con_art5 * 
##     host, family = quasipoisson, data = n1, subset = no_tim > 
##     0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.7302  -1.5775  -0.5168   0.9412   7.4701  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.28015    1.48490  -1.536  0.12864    
## c.mal             0.56452    0.25661   2.200  0.03074 *  
## ln.vol_cub        0.45829    0.08452   5.422 6.22e-07 ***
## cn2              -0.05462    0.04758  -1.148  0.25450    
## c.con_art5        0.07966    0.02794   2.851  0.00556 ** 
## host1             0.80262    0.40694   1.972  0.05207 .  
## c.con_art5:host1 -0.06513    0.04052  -1.607  0.11196    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.261526)
## 
##     Null deviance: 624.76  on 85  degrees of freedom
## Residual deviance: 360.09  on 79  degrees of freedom
##   (8 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an2.f <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.f)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host, 
##     family = quasipoisson, data = n1, subset = no_tim > 0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.6751  -1.6447  -0.7855   0.8684   7.8261  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.07216    1.02646  -2.993   0.0036 ** 
## c.mal             0.49991    0.23689   2.110   0.0377 *  
## ln.vol_cub        0.43164    0.08042   5.367 6.55e-07 ***
## c.con_art5        0.06731    0.02611   2.578   0.0116 *  
## host1             0.47820    0.21732   2.200   0.0304 *  
## c.con_art5:host1 -0.06740    0.03701  -1.821   0.0720 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
## 
##     Null deviance: 652.94  on 92  degrees of freedom
## Residual deviance: 407.37  on 87  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# partial r-sq
an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
           family=quasipoisson, data=n1, subset = host=="C")
summary(an2.1) # 37.6% R-sq
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson, 
##     data = n1, subset = host == "C")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.6270  -1.5586  -0.7073   0.7373   7.9524  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.05242    1.26538  -2.412   0.0189 *  
## c.mal        0.71362    0.34010   2.098   0.0401 *  
## ln.vol_cub   0.42820    0.09911   4.320 5.95e-05 ***
## c.con_art5   0.07210    0.02852   2.528   0.0141 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 6.078791)
## 
##     Null deviance: 512.65  on 63  degrees of freedom
## Residual deviance: 311.59  on 60  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
n1.p <- n1[complete.cases(n1$no_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]
n1.p$no_art_5 <- sqrt(n1.p$no_art_5)
# to drop levels properly, keep vars of interest
n1.p <- n1.p[,c("no_art_5", "c.mal", "ln.vol_cub", "c.con_art5", "host", "no_art_5")]

an.full <- lm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
           data=n1.p, subset = host=="C")
an.art_mal <- lm(no_art_5 ~ ln.vol_cub + c.con_art5, 
           data=n1.p, subset=host=="C")
an.art_vol <- lm(no_art_5 ~ c.mal + c.con_art5, 
           data=n1.p, subset=host=="C")
an.art_con <- lm(no_art_5 ~ c.mal + ln.vol_cub, 
           data=n1.p, subset=host=="C")

an.mal <- lm(c.mal ~ ln.vol_cub + c.con_art5, 
           data=n1.p, subset=host=="C")
an.vol <- lm(ln.vol_cub ~ c.mal + c.con_art5, 
           data=n1.p, subset=host=="C")
an.con <- lm(c.con_art5 ~ c.mal + ln.vol_cub, 
           data=n1.p, subset=host=="C")

r.art_mal <- resid(an.art_mal)
r.art_vol <- resid(an.art_vol)
r.art_con <- resid(an.art_con)
r.mal <- resid(an.mal)
r.vol <- resid(an.vol)
r.con <- resid(an.con)

an.p_mal <- lm(r.art_mal ~ r.mal)
summary(an.p_mal) # r-sq = 5.91 %
## 
## Call:
## lm(formula = r.art_mal ~ r.mal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9844 -0.7116 -0.0408  0.5926  3.6861 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 2.776e-17  1.436e-01   0.000   1.0000  
## r.mal       9.371e-01  4.210e-01   2.226   0.0297 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared:  0.074,  Adjusted R-squared:  0.05906 
## F-statistic: 4.955 on 1 and 62 DF,  p-value: 0.02967
an.p_vol <- lm(r.art_vol ~ r.vol)
summary(an.p_vol) # r-sq = 33.7 %
## 
## Call:
## lm(formula = r.art_vol ~ r.vol)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9844 -0.7116 -0.0408  0.5926  3.6861 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.327e-17  1.436e-01   0.000        1    
## r.vol        5.473e-01  9.519e-02   5.749 2.95e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared:  0.3477, Adjusted R-squared:  0.3372 
## F-statistic: 33.05 on 1 and 62 DF,  p-value: 2.945e-07
an.p_con <- lm(r.art_con ~ r.con)
summary(an.p_con) # r-sq = 7.11 %
## 
## Call:
## lm(formula = r.art_con ~ r.con)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9844 -0.7116 -0.0408  0.5926  3.6861 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -3.174e-18  1.436e-01   0.000   1.0000  
## r.con        1.143e-01  4.739e-02   2.413   0.0188 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared:  0.08584,    Adjusted R-squared:  0.07109 
## F-statistic: 5.822 on 1 and 62 DF,  p-value: 0.0188
# herbivore abundance
anh <- glm(no_herb_5 ~ host*c.con_art5 + cn2 + ln.vol_cub, 
           family=quasipoisson, data=n1)
summary(anh)
## 
## Call:
## glm(formula = no_herb_5 ~ host * c.con_art5 + cn2 + ln.vol_cub, 
##     family = quasipoisson, data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.0557  -1.7191  -0.4770   0.5317   6.5428  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -4.301610   0.919498  -4.678 7.11e-06 ***
## host1             0.499453   0.276353   1.807  0.07301 .  
## c.con_art5        0.050709   0.019107   2.654  0.00894 ** 
## cn2              -0.036775   0.028114  -1.308  0.19315    
## ln.vol_cub        0.609810   0.057103  10.679  < 2e-16 ***
## host1:c.con_art5 -0.008761   0.026901  -0.326  0.74520    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 4.263334)
## 
##     Null deviance: 1314.29  on 136  degrees of freedom
## Residual deviance:  518.29  on 131  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# on arthropod diversity (5 mm)

contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)

## for neg freq dep mal
an3 <- glm(rich_art_5 ~ c.nf_mal_k2.0 * host + ln.vol_cub, 
           family=quasipoisson, data=n1, subset=no_tim > 0)
summary(an3)
## 
## Call:
## glm(formula = rich_art_5 ~ c.nf_mal_k2.0 * host + ln.vol_cub, 
##     family = quasipoisson, data = n1, subset = no_tim > 0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5840  -0.9553  -0.3483   0.7788   2.3267  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -2.72255    0.70002  -3.889 0.000195 ***
## c.nf_mal_k2.0        0.45476    0.19918   2.283 0.024827 *  
## host1               -0.06453    0.16319  -0.395 0.693484    
## ln.vol_cub           0.35269    0.05493   6.421 6.67e-09 ***
## c.nf_mal_k2.0:host1 -0.64648    0.45049  -1.435 0.154820    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.524261)
## 
##     Null deviance: 234.49  on 92  degrees of freedom
## Residual deviance: 140.23  on 88  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
##


an3 <- glm(rich_art_5 ~ c.mal*host + ln.vol_cub + host * c.con_art5, 
           family=quasipoisson, data=n1, subset=no_tim > 0)
summary(an3)
## 
## Call:
## glm(formula = rich_art_5 ~ c.mal * host + ln.vol_cub + host * 
##     c.con_art5, family = quasipoisson, data = n1, subset = no_tim > 
##     0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3865  -1.0153  -0.2121   0.7193   2.0200  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.75709    0.68817  -4.006 0.000131 ***
## c.mal             0.72859    0.21402   3.404 0.001010 ** 
## host1            -0.02279    0.17607  -0.129 0.897303    
## ln.vol_cub        0.35623    0.05413   6.581 3.51e-09 ***
## c.con_art5        0.05687    0.01870   3.041 0.003126 ** 
## c.mal:host1      -0.98839    0.39786  -2.484 0.014921 *  
## host1:c.con_art5 -0.08287    0.02862  -2.896 0.004791 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
## 
##     Null deviance: 234.49  on 92  degrees of freedom
## Residual deviance: 122.51  on 86  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an3 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + host*c.con_art5, 
           family=quasipoisson, data=n1, subset=no_tim_gs > 0)
summary(an3)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub + 
##     host * c.con_art5, family = quasipoisson, data = n1, subset = no_tim_gs > 
##     0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3865  -1.0153  -0.2121   0.7193   2.0200  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.75709    0.68817  -4.006 0.000131 ***
## host1            -0.02279    0.17607  -0.129 0.897303    
## c.mal             0.72859    0.21402   3.404 0.001010 ** 
## ln.vol_cub        0.35623    0.05413   6.581 3.51e-09 ***
## c.con_art5        0.05687    0.01870   3.041 0.003126 ** 
## host1:c.mal      -0.98839    0.39786  -2.484 0.014921 *  
## host1:c.con_art5 -0.08287    0.02862  -2.896 0.004791 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
## 
##     Null deviance: 234.49  on 92  degrees of freedom
## Residual deviance: 122.51  on 86  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
contrasts(n1$host) <- c(0,1)
an3.1 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + c.con_art5, 
             family=quasipoisson, data=n1, subset=no_tim>1)
summary(an3.1)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub + 
##     c.con_art5, family = quasipoisson, data = n1, subset = no_tim > 
##     1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5165  -1.1008  -0.1576   0.9107   2.0497  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.50369    0.82816  -3.023  0.00368 ** 
## host1        0.61763    0.25225   2.448  0.01729 *  
## c.mal       -0.83967    0.61093  -1.374  0.17442    
## ln.vol_cub   0.28991    0.06540   4.433 4.02e-05 ***
## c.con_art5   0.04197    0.01804   2.326  0.02341 *  
## host1:c.mal  1.54806    0.67284   2.301  0.02490 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.446096)
## 
##     Null deviance: 153.997  on 65  degrees of freedom
## Residual deviance:  90.206  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an3.2 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + host*c.con_art5 +
               cn_rat, 
           family=quasipoisson, data=n1, subset=no_tim_gs > 3)
summary(an3.2)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub + 
##     host * c.con_art5 + cn_rat, family = quasipoisson, data = n1, 
##     subset = no_tim_gs > 3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6129  -0.9568  -0.2084   0.6939   2.2196  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      -0.655505   2.317628  -0.283   0.7791  
## host1            -0.390998   0.604090  -0.647   0.5221  
## c.mal            -0.093110   1.124048  -0.083   0.9345  
## ln.vol_cub        0.327657   0.124026   2.642   0.0127 *
## c.con_art5        0.003396   0.053773   0.063   0.9500  
## cn_rat           -0.066657   0.056740  -1.175   0.2487  
## host1:c.mal       1.325556   1.247692   1.062   0.2960  
## host1:c.con_art5  0.071246   0.062926   1.132   0.2660  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.635073)
## 
##     Null deviance: 84.078  on 39  degrees of freedom
## Residual deviance: 52.518  on 32  degrees of freedom
##   (2 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# final model

an3.1 <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + cn2 + host*c.con_art5, 
           family=quasipoisson, data=n1)
summary(an3.1)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + cn2 + 
##     host * c.con_art5, family = quasipoisson, data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3927  -1.0121  -0.2293   0.6199   2.1031  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.17357    1.18643  -1.832  0.07076 .  
## host1            -0.19521    0.31487  -0.620  0.53708    
## c.mal            -0.30162    0.34848  -0.866  0.38940    
## ln.vol_cub        0.37460    0.05774   6.487 7.32e-09 ***
## cn2              -0.02960    0.03520  -0.841  0.40301    
## c.con_art5       -0.02420    0.02222  -1.089  0.27959    
## host1:c.mal       1.10013    0.40981   2.684  0.00887 ** 
## host1:c.con_art5  0.09438    0.03096   3.049  0.00314 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.346103)
## 
##     Null deviance: 221.89  on 85  degrees of freedom
## Residual deviance: 110.07  on 78  degrees of freedom
##   (60 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an3.f <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + host*c.con_art5, 
           family=quasipoisson, data=n1)
summary(an3.f)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + host * 
##     c.con_art5, family = quasipoisson, data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3865  -1.0153  -0.2121   0.7193   2.0200  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.77989    0.63926  -4.349 3.75e-05 ***
## host1             0.02279    0.17607   0.129  0.89730    
## c.mal            -0.25980    0.33514  -0.775  0.44034    
## ln.vol_cub        0.35623    0.05413   6.581 3.51e-09 ***
## c.con_art5       -0.02600    0.02146  -1.212  0.22897    
## host1:c.mal       0.98839    0.39786   2.484  0.01492 *  
## host1:c.con_art5  0.08287    0.02862   2.896  0.00479 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
## 
##     Null deviance: 234.49  on 92  degrees of freedom
## Residual deviance: 122.51  on 86  degrees of freedom
##   (53 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
# partial r-sq

an3 <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, 
           family=quasipoisson, data=n1, subset=host=="C")
summary(an3) # R-sq = 53% total
## 
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson, 
##     data = n1, subset = host == "C")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3858  -1.1404  -0.2420   0.8255   2.0197  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.75955    0.79599  -3.467  0.00098 ***
## c.mal        0.72862    0.22775   3.199  0.00220 ** 
## ln.vol_cub   0.35643    0.06266   5.688 4.05e-07 ***
## c.con_art5   0.05688    0.01992   2.855  0.00590 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.538553)
## 
##     Null deviance: 184.732  on 63  degrees of freedom
## Residual deviance:  97.756  on 60  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
n1.p <- n1[complete.cases(n1$rich_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]
n1.p$rich_art_5 <- sqrt(n1.p$rich_art_5)

an.rich_mal <- lm(rich_art_5 ~ ln.vol_cub + c.con_art5, 
           data=n1.p, subset=host=="C")
an.rich_vol <- lm(rich_art_5 ~ c.mal + c.con_art5, 
           data=n1.p, subset=host=="C")
an.rich_con <- lm(rich_art_5 ~ c.mal + ln.vol_cub, 
           data=n1.p, subset=host=="C")
an.mal <- lm(c.mal ~ ln.vol_cub + c.con_art5, 
           data=n1.p, subset=host=="C")
an.vol <- lm(ln.vol_cub ~ c.mal + c.con_art5, 
           data=n1.p, subset=host=="C")
an.con <- lm(c.con_art5 ~ c.mal + ln.vol_cub, 
           data=n1.p, subset=host=="C")

r.rich_mal <- resid(an.rich_mal)
r.rich_vol <- resid(an.rich_vol)
r.rich_con <- resid(an.rich_con)
r.mal <- resid(an.mal)
r.vol <- resid(an.vol)
r.con <- resid(an.con)

an.p_mal <- lm(r.rich_mal ~ r.mal)
summary(an.p_mal) # r-sq = 9.92 %
## 
## Call:
## lm(formula = r.rich_mal ~ r.mal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.55867 -0.48803  0.00644  0.56164  1.30531 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -3.816e-17  8.486e-02   0.000  1.00000   
## r.mal        7.012e-01  2.488e-01   2.818  0.00648 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared:  0.1135, Adjusted R-squared:  0.09924 
## F-statistic: 7.941 on 1 and 62 DF,  p-value: 0.006476
an.p_vol <- lm(r.rich_vol ~ r.vol)
summary(an.p_vol) # r-sq = 42.3 %
## 
## Call:
## lm(formula = r.rich_vol ~ r.vol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.55867 -0.48803  0.00644  0.56164  1.30531 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.469e-17  8.486e-02   0.000        1    
## r.vol       3.864e-01  5.626e-02   6.868 3.63e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared:  0.4321, Adjusted R-squared:  0.4229 
## F-statistic: 47.18 on 1 and 62 DF,  p-value: 3.631e-09
an.p_con <- lm(r.rich_con ~ r.con)
summary(an.p_con) # r-sq = 5.84 %
## 
## Call:
## lm(formula = r.rich_con ~ r.con)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.55867 -0.48803  0.00644  0.56164  1.30531 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1.907e-17  8.486e-02   0.000   1.0000  
## r.con        6.205e-02  2.801e-02   2.215   0.0304 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared:  0.07335,    Adjusted R-squared:  0.0584 
## F-statistic: 4.908 on 1 and 62 DF,  p-value: 0.03042
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)

an5 <- lm(cn2 ~ c.mal + host + ln.vol_cub , data=n1, subset=no_tim > 0)
summary(an5)
## 
## Call:
## lm(formula = cn2 ~ c.mal + host + ln.vol_cub, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3563 -1.3177  0.0028  1.1640  4.5875 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 19.94660    1.90134  10.491   <2e-16 ***
## c.mal       -1.31193    0.54549  -2.405   0.0184 *  
## host1        7.56694    0.46137  16.401   <2e-16 ***
## ln.vol_cub   0.05646    0.15487   0.365   0.7164    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.81 on 82 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.8068, Adjusted R-squared:  0.7997 
## F-statistic: 114.1 on 3 and 82 DF,  p-value: < 2.2e-16
an6 <- lm(cn2 ~ pC, data=n1)
an7 <- lm(cn2 ~ pN, data=n1)
summary(an6)
## 
## Call:
## lm(formula = cn2 ~ pC, data = n1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.898 -3.701 -1.085  3.810 14.042 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  27.0853     3.2998   8.208 1.59e-13 ***
## pC           -0.2731     0.3185  -0.858    0.393    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.656 on 135 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.005417,   Adjusted R-squared:  -0.00195 
## F-statistic: 0.7353 on 1 and 135 DF,  p-value: 0.3927
summary(an7)
## 
## Call:
## lm(formula = cn2 ~ pN, data = n1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0101 -1.7586 -0.3754  1.3619 18.4994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   39.854      1.030   38.70   <2e-16 ***
## pN           -35.430      2.279  -15.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.795 on 135 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.6417, Adjusted R-squared:  0.639 
## F-statistic: 241.7 on 1 and 135 DF,  p-value: < 2.2e-16
contrasts(n1$host) <- c(0, 1)
an5 <- lm(pC ~ c.mal + host, data=n1, subset=no_tim>0)
summary(an5)
## 
## Call:
## lm(formula = pC ~ c.mal + host, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0288 -0.5477 -0.0229  0.5032 11.0203 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.1707     0.2980  34.132   <2e-16 ***
## c.mal         0.2730     0.4463   0.612    0.542    
## host1         0.2438     0.3606   0.676    0.501    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.482 on 83 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.01341,    Adjusted R-squared:  -0.01036 
## F-statistic: 0.564 on 2 and 83 DF,  p-value: 0.5711
contrasts(n1$host) <- c(0, 1)
an5 <- lm(pN ~ c.mal + host + no_herb , data=n1, subset=no_tim>0)
summary(an5)
## 
## Call:
## lm(formula = pN ~ c.mal + host + no_herb, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15346 -0.03789  0.00174  0.02453  0.58215 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.3535096  0.0202036  17.497  < 2e-16 ***
## c.mal       0.0380642  0.0251261   1.515    0.134    
## host1       0.1418512  0.0199923   7.095 4.19e-10 ***
## no_herb     0.0008782  0.0008526   1.030    0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08202 on 82 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.4475, Adjusted R-squared:  0.4273 
## F-statistic: 22.14 on 3 and 82 DF,  p-value: 1.35e-10
# final model
an5.1 <- lm(cn2 ~ c.mal + host + ln.vol_cub + no_art_5, data=n1, subset=no_tim > 0)
summary(an5.1)
## 
## Call:
## lm(formula = cn2 ~ c.mal + host + ln.vol_cub + no_art_5, data = n1, 
##     subset = no_tim > 0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3052 -1.1541 -0.0549  1.0570  4.3339 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.37517    1.91214  13.794   <2e-16 ***
## c.mal       -1.13413    0.55220  -2.054   0.0432 *  
## host1       -7.71009    0.46613 -16.541   <2e-16 ***
## ln.vol_cub   0.19395    0.17646   1.099   0.2750    
## no_art_5    -0.04080    0.02585  -1.579   0.1183    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.794 on 81 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.8125, Adjusted R-squared:  0.8033 
## F-statistic: 87.76 on 4 and 81 DF,  p-value: < 2.2e-16
an5.2 <- lm(cn2 ~ c.mal + host + no_art_5, data=n1, subset=no_tim > 0)
summary(an5.2)# R-sq = 80.2 %
## 
## Call:
## lm(formula = cn2 ~ c.mal + host + no_art_5, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.377 -1.086 -0.082  1.093  4.516 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 28.42373    0.42761  66.471   <2e-16 ***
## c.mal       -1.17752    0.55149  -2.135   0.0357 *  
## host1       -7.53107    0.43730 -17.222   <2e-16 ***
## no_art_5    -0.02678    0.02251  -1.190   0.2375    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.796 on 82 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.8097, Adjusted R-squared:  0.8028 
## F-statistic: 116.3 on 3 and 82 DF,  p-value: < 2.2e-16
an5.f <- lm(cn2 ~ c.mal + host, data=n1, subset=no_tim > 0)
summary(an5.f) 
## 
## Call:
## lm(formula = cn2 ~ c.mal + host, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.377 -1.280 -0.009  1.170  4.624 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.1514     0.3621  77.739   <2e-16 ***
## c.mal        -1.3051     0.5423  -2.407   0.0183 *  
## host1        -7.5170     0.4382 -17.153   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.8 on 83 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.8064, Adjusted R-squared:  0.8018 
## F-statistic: 172.9 on 2 and 83 DF,  p-value: < 2.2e-16
## neg freq dep mal
an5.f <- lm(cn2 ~ c.nf_mal_k0.5 + host, data=n1, subset=no_tim > 0)
summary(an5.f)
## 
## Call:
## lm(formula = cn2 ~ c.nf_mal_k0.5 + host, data = n1, subset = no_tim > 
##     0)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1498 -1.3665 -0.0076  1.2264  4.2070 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    28.2882     0.3566  79.330   <2e-16 ***
## c.nf_mal_k0.5  -1.1872     0.5567  -2.133   0.0359 *  
## host1          -7.7235     0.4270 -18.086   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.813 on 83 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.8037, Adjusted R-squared:  0.799 
## F-statistic: 169.9 on 2 and 83 DF,  p-value: < 2.2e-16
##

# partial r-sq
p.n1 <- n1[complete.cases(n1$host, n1$c.mal, n1$cn2), ]
p.n1$c.host <- ifelse(p.n1$host == "A", -.5, .5)

an.mal <- lm(c.mal ~ c.host, data=p.n1)
an.host <- lm(c.host ~ c.mal, data=p.n1)
an.cn_h <- lm(cn2 ~ c.host, data=p.n1)
an.cn_m <- lm(cn2 ~ c.mal, data=p.n1)

r.mal <- resid(an.mal)
r.host <- resid(an.host)
r.cn_h <- resid(an.cn_h)
r.cn_m <- resid(an.cn_m)

an.p_mal <- lm(r.cn_h ~ r.mal) 
summary(an.p_mal) # r-sq = 5.41%
## 
## Call:
## lm(formula = r.cn_h ~ r.mal)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.377 -1.280 -0.009  1.170  4.624 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -4.202e-17  1.930e-01   0.000   1.0000  
## r.mal       -1.305e+00  5.391e-01  -2.421   0.0176 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.79 on 84 degrees of freedom
## Multiple R-squared:  0.06523,    Adjusted R-squared:  0.0541 
## F-statistic: 5.861 on 1 and 84 DF,  p-value: 0.01763
an.p_host <- lm(r.cn_m ~ r.host)
summary(an.p_host) # r-sq = 77.7%
## 
## Call:
## lm(formula = r.cn_m ~ r.host)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.377 -1.280 -0.009  1.170  4.624 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.551e-16  1.930e-01    0.00        1    
## r.host      -7.517e+00  4.356e-01  -17.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.79 on 84 degrees of freedom
## Multiple R-squared:   0.78,  Adjusted R-squared:  0.7773 
## F-statistic: 297.8 on 1 and 84 DF,  p-value: < 2.2e-16
# on biomass
# changed starts from 0
an2 <- glm(biomass ~ mal, family=Gamma, start=c(0.0001,0.0001),
           data=n1[n1$biomass > 0,])
summary(an2)
## 
## Call:
## glm(formula = biomass ~ mal, family = Gamma, data = n1[n1$biomass > 
##     0, ], start = c(1e-04, 1e-04))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4097  -1.0986  -0.4272   0.2647   2.0929  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.295      1.240   5.885 7.47e-08 ***
## mal           -3.431      1.733  -1.980   0.0509 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 1.196149)
## 
##     Null deviance: 111.56  on 87  degrees of freedom
## Residual deviance: 106.84  on 86  degrees of freedom
##   (41 observations deleted due to missingness)
## AIC: -120.38
## 
## Number of Fisher Scoring iterations: 21
## connectivity on maladaptation?

mal_nomal <- cbind(n1_occ$no_mal, n1_occ$no_adp)
grn_str <- cbind(n1_occ$no_grn_tim, n1_occ$no_str_tim)

#an1 <- glm(tims ~ n1_occ$con_tim, family=binomial) # no effect of connectivity, !!what is tims
an2 <- glm(mal ~ con_tim, data=n1_occ)
summary(an2)
## 
## Call:
## glm(formula = mal ~ con_tim, data = n1_occ)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.49278  -0.43417   0.02203   0.29194   0.56315  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.49289    0.05898   8.357 6.96e-13 ***
## con_tim     -0.01668    0.02085  -0.800    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1430296)
## 
##     Null deviance: 13.107  on 92  degrees of freedom
## Residual deviance: 13.016  on 91  degrees of freedom
## AIC: 87.043
## 
## Number of Fisher Scoring iterations: 2
an3 <- glm(mal_nomal ~ no_art_5, family=binomial, data=n1_occ, subset=no_tim>0)

summary(an3)
## 
## Call:
## glm(formula = mal_nomal ~ no_art_5, family = binomial, data = n1_occ, 
##     subset = no_tim > 0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5271  -1.0306  -0.3327   0.8986   2.5854  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.338683   0.166728  -2.031   0.0422 *
## no_art_5    -0.002842   0.010716  -0.265   0.7909  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 133.15  on 92  degrees of freedom
## Residual deviance: 133.08  on 91  degrees of freedom
## AIC: 239.92
## 
## Number of Fisher Scoring iterations: 3
an4 <- lm(mal ~ no_art_5 + ln.vol_cub + host + c.con_art5, data=n1)
summary(an4)
## 
## Call:
## lm(formula = mal ~ no_art_5 + ln.vol_cub + host + c.con_art5, 
##     data = n1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55503 -0.27132 -0.03773  0.27866  0.76125 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.515242   0.330751   1.558  0.12287   
## no_art_5     0.009593   0.004962   1.933  0.05643 . 
## ln.vol_cub  -0.027387   0.030703  -0.892  0.37482   
## host1        0.248194   0.094221   2.634  0.00996 **
## c.con_art5  -0.008293   0.009454  -0.877  0.38281   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3568 on 88 degrees of freedom
##   (53 observations deleted due to missingness)
## Multiple R-squared:  0.1451, Adjusted R-squared:  0.1062 
## F-statistic: 3.734 on 4 and 88 DF,  p-value: 0.007463
# plant growth

hist(n1$avg_growth)

an1 <- lm(avg_growth ~  con_art5 + host + no_art_5, data=n1)
summary(an1)
## 
## Call:
## lm(formula = avg_growth ~ con_art5 + host + no_art_5, data = n1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.527 -12.393  -1.892  10.114  45.407 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  81.8856     4.0065  20.438   <2e-16 ***
## con_art5     -0.9006     0.3759  -2.396   0.0183 *  
## host1       -40.4208     3.9231 -10.303   <2e-16 ***
## no_art_5     -0.5023     0.2250  -2.232   0.0277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.14 on 105 degrees of freedom
##   (37 observations deleted due to missingness)
## Multiple R-squared:  0.5578, Adjusted R-squared:  0.5452 
## F-statistic: 44.15 on 3 and 105 DF,  p-value: < 2.2e-16
hist(n1$avg_no_flwrs)

n1$no_not_herb_5 <- n1$no_art_5 - n1$no_herb_5
an2 <- glm(avg_no_flwrs ~ no_herb_5 + ln.vol_cub, data=n1, family=quasipoisson)
summary(an2)
## 
## Call:
## glm(formula = avg_no_flwrs ~ no_herb_5 + ln.vol_cub, family = quasipoisson, 
##     data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8752  -1.0093  -0.5487   0.2925   2.5442  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -3.53768    1.22170  -2.896  0.00508 **
## no_herb_5    0.01450    0.01187   1.221  0.22618   
## ln.vol_cub   0.29910    0.10407   2.874  0.00540 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.286064)
## 
##     Null deviance: 109.804  on 70  degrees of freedom
## Residual deviance:  83.553  on 68  degrees of freedom
##   (75 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)

an.s1 <- lm(mas ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim + c.con_art5, data=n1[(n1$no_art_5 > 2 & n1$no_tim > 0 & n1$c.mal < 0.6), ])
summary(an.s1)
## 
## Call:
## lm(formula = mas ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim + c.con_art5, data = n1[(n1$no_art_5 > 2 & n1$no_tim > 
##     0 & n1$c.mal < 0.6), ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60977 -0.18461 -0.01867  0.15665  0.64327 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.681677   0.053586 -12.721  < 2e-16 ***
## host1             0.214007   0.103443   2.069  0.04242 *  
## c.mal             0.150830   0.129742   1.163  0.24914    
## c.no_art_5       -0.030005   0.004684  -6.406 1.72e-08 ***
## c.no_tim         -0.014150   0.009853  -1.436  0.15560    
## c.con_art5        0.027374   0.008206   3.336  0.00139 ** 
## host1:c.mal       0.116187   0.228700   0.508  0.61310    
## host1:c.no_art_5 -0.021955   0.010567  -2.078  0.04157 *  
## c.mal:c.no_tim    0.101103   0.036701   2.755  0.00755 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2837 on 67 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5893, Adjusted R-squared:  0.5403 
## F-statistic: 12.02 on 8 and 67 DF,  p-value: 1.833e-10
an.sf <- lm(mas ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim + c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an.sf)
## 
## Call:
## lm(formula = mas ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim + c.con_art5, data = n1[n1$no_art_5 > 2 & n1$no_tim > 
##     0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60977 -0.18461 -0.01867  0.15665  0.64327 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.681677   0.053586 -12.721  < 2e-16 ***
## host1             0.214007   0.103443   2.069  0.04242 *  
## c.mal             0.150830   0.129742   1.163  0.24914    
## c.no_art_5       -0.030005   0.004684  -6.406 1.72e-08 ***
## c.no_tim         -0.014150   0.009853  -1.436  0.15560    
## c.con_art5        0.027374   0.008206   3.336  0.00139 ** 
## host1:c.mal       0.116187   0.228700   0.508  0.61310    
## host1:c.no_art_5 -0.021955   0.010567  -2.078  0.04157 *  
## c.mal:c.no_tim    0.101103   0.036701   2.755  0.00755 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2837 on 67 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5893, Adjusted R-squared:  0.5403 
## F-statistic: 12.02 on 8 and 67 DF,  p-value: 1.833e-10
## neg freq dep mal
#an.sf <- lm(mas ~ host*c.nf_mal_k2.0 + host*c.no_art_5 + c.nf_mal_k2.0*c.no_tim + c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
#summary(an.sf)
###

an1 <- lm(mai1 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an1)
## 
## Call:
## lm(formula = mai1 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.50939 -0.40094 -0.02639  0.31742  2.32790 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.488093   0.115158   4.238 6.94e-05 ***
## host1            -0.507877   0.212174  -2.394  0.01945 *  
## c.mal             0.054880   0.283720   0.193  0.84720    
## c.no_art_5        0.052519   0.010370   5.065 3.34e-06 ***
## c.no_tim          0.004581   0.021266   0.215  0.83010    
## host1:c.mal      -1.054458   0.502242  -2.100  0.03949 *  
## host1:c.no_art_5  0.044190   0.023411   1.888  0.06335 .  
## c.mal:c.no_tim   -0.276703   0.080579  -3.434  0.00102 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6286 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4525, Adjusted R-squared:  0.3961 
## F-statistic: 8.028 on 7 and 68 DF,  p-value: 4.402e-07
an20 <- lm(mai20 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an20)
## 
## Call:
## lm(formula = mai20 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3788 -0.4707 -0.0875  0.3971  3.4183 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.06086    0.17739  -0.343 0.732602    
## host1             1.22800    0.32684   3.757 0.000359 ***
## c.mal             0.03294    0.43705   0.075 0.940145    
## c.no_art_5       -0.01297    0.01597  -0.812 0.419549    
## c.no_tim          0.02004    0.03276   0.612 0.542669    
## host1:c.mal       2.45150    0.77366   3.169 0.002295 ** 
## host1:c.no_art_5 -0.06586    0.03606  -1.826 0.072183 .  
## c.mal:c.no_tim    0.33483    0.12413   2.698 0.008801 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9683 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2636, Adjusted R-squared:  0.1878 
## F-statistic: 3.477 on 7 and 68 DF,  p-value: 0.003047
an3 <- lm(mai3 ~ host*c.mal+ host*c.no_art_5 + c.mal*c.no_tim, 
          data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an3)
## 
## Call:
## lm(formula = mai3 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99069 -0.30888 -0.04333  0.22945  1.86612 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.430309   0.088111   4.884 6.62e-06 ***
## host1            -0.325153   0.162341  -2.003 0.049177 *  
## c.mal             0.052570   0.217083   0.242 0.809381    
## c.no_art_5        0.045625   0.007934   5.750 2.31e-07 ***
## c.no_tim          0.006208   0.016271   0.382 0.703982    
## host1:c.mal      -0.685410   0.384281  -1.784 0.078948 .  
## host1:c.no_art_5  0.032606   0.017913   1.820 0.073120 .  
## c.mal:c.no_tim   -0.212331   0.061654  -3.444 0.000986 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4809 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4951, Adjusted R-squared:  0.4432 
## F-statistic: 9.528 on 7 and 68 DF,  p-value: 3.449e-08
an5 <- lm(mai5 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
          data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an5)
## 
## Call:
## lm(formula = mai5 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5089 -0.2210 -0.0606  0.1439  1.4043 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.372524   0.063327   5.883 1.36e-07 ***
## host1            -0.142429   0.116677  -1.221  0.22641    
## c.mal             0.050260   0.156021   0.322  0.74834    
## c.no_art_5        0.038731   0.005703   6.792 3.36e-09 ***
## c.no_tim          0.007836   0.011694   0.670  0.50508    
## host1:c.mal      -0.316361   0.276189  -1.145  0.25604    
## host1:c.no_art_5  0.021021   0.012874   1.633  0.10713    
## c.mal:c.no_tim   -0.147959   0.044311  -3.339  0.00137 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3457 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.558,  Adjusted R-squared:  0.5125 
## F-statistic: 12.27 on 7 and 68 DF,  p-value: 4.965e-10
an7 <- lm(mai7 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
          data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an7)
## 
## Call:
## lm(formula = mai7 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38411 -0.15051 -0.05102  0.08080  0.94254 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.314740   0.044738   7.035 1.23e-09 ***
## host1             0.040294   0.082429   0.489  0.62653    
## c.mal             0.047951   0.110224   0.435  0.66492    
## c.no_art_5        0.031837   0.004029   7.903 3.29e-11 ***
## c.no_tim          0.009464   0.008262   1.145  0.25602    
## host1:c.mal       0.052687   0.195119   0.270  0.78796    
## host1:c.no_art_5  0.009436   0.009095   1.037  0.30318    
## c.mal:c.no_tim   -0.083587   0.031305  -2.670  0.00948 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2442 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.623,  Adjusted R-squared:  0.5842 
## F-statistic: 16.05 on 7 and 68 DF,  p-value: 2.902e-12
an8 <- lm(mai8 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
          data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an8)
## 
## Call:
## lm(formula = mai8 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33478 -0.16438 -0.03142  0.10841  0.74215 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.285848   0.040637   7.034 1.23e-09 ***
## host1             0.131656   0.074872   1.758   0.0832 .  
## c.mal             0.046796   0.100119   0.467   0.6417    
## c.no_art_5        0.028390   0.003659   7.758 6.02e-11 ***
## c.no_tim          0.010277   0.007504   1.370   0.1753    
## host1:c.mal       0.237211   0.177231   1.338   0.1852    
## host1:c.no_art_5  0.003644   0.008261   0.441   0.6606    
## c.mal:c.no_tim   -0.051401   0.028435  -1.808   0.0751 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2218 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.6229, Adjusted R-squared:  0.584 
## F-statistic: 16.04 on 7 and 68 DF,  p-value: 2.933e-12
an9 <- lm(mai9 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
          data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an9)
## 
## Call:
## lm(formula = mai9 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38322 -0.15543 -0.04908  0.11741  0.70110 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.256956   0.041681   6.165 4.37e-08 ***
## host1             0.223018   0.076796   2.904  0.00496 ** 
## c.mal             0.045641   0.102692   0.444  0.65813    
## c.no_art_5        0.024943   0.003753   6.646 6.13e-09 ***
## c.no_tim          0.011091   0.007697   1.441  0.15418    
## host1:c.mal       0.421735   0.181785   2.320  0.02335 *  
## host1:c.no_art_5 -0.002149   0.008474  -0.254  0.80060    
## c.mal:c.no_tim   -0.019215   0.029165  -0.659  0.51222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2275 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.5745, Adjusted R-squared:  0.5307 
## F-statistic: 13.12 on 7 and 68 DF,  p-value: 1.461e-10
an10 <- lm(mai10 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
           data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an10)
## 
## Call:
## lm(formula = mai10 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45265 -0.16372 -0.05317  0.14843  0.82478 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.228063   0.047533   4.798 9.13e-06 ***
## host1             0.314380   0.087578   3.590  0.00062 ***
## c.mal             0.044486   0.117110   0.380  0.70523    
## c.no_art_5        0.021497   0.004280   5.022 3.93e-06 ***
## c.no_tim          0.011905   0.008778   1.356  0.17950    
## host1:c.mal       0.606259   0.207308   2.924  0.00468 ** 
## host1:c.no_art_5 -0.007941   0.009663  -0.822  0.41409    
## c.mal:c.no_tim    0.012971   0.033260   0.390  0.69778    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2595 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4891, Adjusted R-squared:  0.4365 
## F-statistic:   9.3 on 7 and 68 DF,  p-value: 5.018e-08
an18 <- lm(mai18 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
           data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an18)
## 
## Call:
## lm(formula = mai18 ~ host * c.mal + host * c.no_art_5 + c.mal * 
##     c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.98226 -0.39681 -0.08116  0.33437  2.89961 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.003074   0.148785  -0.021 0.983578    
## host1             1.045274   0.274130   3.813 0.000298 ***
## c.mal             0.035248   0.366569   0.096 0.923679    
## c.no_art_5       -0.006079   0.013398  -0.454 0.651466    
## c.no_tim          0.018416   0.027476   0.670 0.504960    
## host1:c.mal       2.082452   0.648900   3.209 0.002032 ** 
## host1:c.no_art_5 -0.054280   0.030247  -1.795 0.077173 .  
## c.mal:c.no_tim    0.270458   0.104109   2.598 0.011491 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8121 on 68 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2619, Adjusted R-squared:  0.1859 
## F-statistic: 3.446 on 7 and 68 DF,  p-value: 0.00325
# predictions on Adenostoma
pred.slope <- predict(an.sf, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE) # changed an.s for an.sf
pred.int <- predict(an1, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)
pred.int20 <- predict(an20, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)


# predictions on Ceanothus - do not plot as not significant
# pred.slope <- predict(an.sf, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE) # changed an.s for an.s1
# pred.int <- predict(an1, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)
# pred.int20 <- predict(an20, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0), c.con_art5=0, se.fit=TRUE)

# partial r-sq
an.s2 <- lm(mas ~ c.mal + c.no_art_5, data=n1[(n1$no_art_5 > 2 & n1$no_tim > 0 & n1$host == "A"), ])
summary(an.s2) # R-sq = 18.4 %
## 
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1[(n1$no_art_5 > 
##     2 & n1$no_tim > 0 & n1$host == "A"), ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62088 -0.20001  0.03397  0.15980  0.85752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.40958    0.10117  -4.048 0.000536 ***
## c.mal        0.18082    0.21821   0.829 0.416207    
## c.no_art_5  -0.04825    0.01169  -4.127 0.000443 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3531 on 22 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.4409, Adjusted R-squared:  0.3901 
## F-statistic: 8.675 on 2 and 22 DF,  p-value: 0.001668
n1.p <- n1[complete.cases(n1$mas, n1$c.mal,n1$c.no_art_5), ]
# to drop levels properly, keep vars of interest
n1.p <- n1.p[,c("c.no_art_5", "c.mal", "no_tim", "host", "no_art_5", "mas")]

an.mas.art <- lm(mas ~ c.no_art_5, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.mas.mal <- lm(mas ~ c.mal, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.mal.art <- lm(c.mal ~ c.no_art_5, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.art.mal <- lm(c.no_art_5 ~ c.mal, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])

r.mas <- resid(an.mas.art)
r.mal <- resid(an.mal.art)
r.art <- resid(an.art.mal)

an.p_mal <- lm(r.mas ~ r.mal)
summary(an.p_mal) # r-sq = 0.4% on Adenostoma
## 
## Call:
## lm(formula = r.mas ~ r.mal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62088 -0.20001  0.03397  0.15980  0.85752 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.353e-18  6.906e-02   0.000    1.000
## r.mal       1.808e-01  2.134e-01   0.847    0.406
## 
## Residual standard error: 0.3453 on 23 degrees of freedom
## Multiple R-squared:  0.03027,    Adjusted R-squared:  -0.0119 
## F-statistic: 0.7179 on 1 and 23 DF,  p-value: 0.4056
an.p_art <- lm(r.mal~ r.art)
summary(an.p_art) # r-sq = 7% on Adenostoma
## 
## Call:
## lm(formula = r.mal ~ r.art)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2155 -0.2155 -0.2155  0.1664  0.6756 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.051e-17  6.370e-02   0.000    1.000
## r.art       -1.768e-02  1.055e-02  -1.676    0.107
## 
## Residual standard error: 0.3185 on 23 degrees of freedom
## Multiple R-squared:  0.1089, Adjusted R-squared:  0.07012 
## F-statistic:  2.81 on 1 and 23 DF,  p-value: 0.1072
pdf(file="figures/mas_curves_by_mal4A.pdf")

plot(x=c(0, 19), y=c(-1, 3), type="n", las=1,
     xlab="rank log biomass", ylab="log abundance", axes=FALSE)
axis(1, at=c(0, 4, 9, 14, 19), labels=c(1, 5, 10, 15, 20))
axis(2, at=seq(-1, 3, by=.5), 
     labels=c("-1", "", "0", "", "1", "", "2", "", "3"), las=1)

# lines(x=c(0, 19), y=c(pred.int$fit[1]-pred.int$se.fit[1], 
#                       pred.int20$fit[1]+pred.int$se.fit[1]),
#       col="blue", lty=2)
# lines(x=c(0, 19), y=c(pred.int$fit[1]+pred.int$se.fit[1], 
#                       pred.int20$fit[1]-pred.int$se.fit[1]),
#       col="blue", lty=2)

polygon(x=c(0, 0, 19, 19), y=c(pred.int$fit[1] - pred.int$se.fit[1], 
                               pred.int$fit[1] + pred.int$se.fit[1], 
                               pred.int20$fit[1] - pred.int20$se.fit[1], 
                               pred.int20$fit[1] + pred.int20$se.fit[1]), 
        lwd=.5, col="blue", density=150)

lines(x=c(0, 19), y=c(pred.int$fit[1], pred.int20$fit[1]),
      col="blue", lwd=2)

# lines(x=c(0, 19), y=c(pred.int$fit[2]-pred.int$se.fit[2], 
#                       pred.int20$fit[2]+pred.int$se.fit[2]),
#       col="red", lty=2)
# lines(x=c(0, 19), y=c(pred.int$fit[2]+pred.int$se.fit[2], 
#                       pred.int20$fit[2]-pred.int$se.fit[2]),
#       col="red", lty=2)

polygon(x=c(0, 0, 19, 19), y=c(pred.int$fit[2] - pred.int$se.fit[2], 
                               pred.int$fit[2] + pred.int$se.fit[2], 
                               pred.int20$fit[2] - pred.int20$se.fit[2], 
                               pred.int20$fit[2] + pred.int20$se.fit[2]), 
        lwd=.5, col="red", density=150)

lines(x=c(0, 19), y=c(pred.int$fit[2], pred.int20$fit[2]),
      col="red", lwd=2)

abline(v=c(2, 8), lty=2, lwd=1, col="darkgreen")

legend(3, 2.5, legend=c("maladapted", "well adapted"), fill=c("red", "blue"), 
       bty="n", cex=.9)

arrows(x0=c(2, 8), x1=c(0, 11), y0=c(0.075, 0.5), y1=c(0.075, 0.5), length=.1)
text(x=c(1, 9.25), y=c(0.2, .6), labels=c("p < 0.05", "p < 0.05"), cex=0.85)

dev.off()
## quartz_off_screen 
##                 2
contrasts(n1$host) <- c(0, 1)

an1 <- lm(mas ~ no_art_5 + c.mal*host, data=n1[n1$mas<0.6 & n1$no_art_5 > 5, ])
summary(an1)
## 
## Call:
## lm(formula = mas ~ no_art_5 + c.mal * host, data = n1[n1$mas < 
##     0.6 & n1$no_art_5 > 5, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54968 -0.17396 -0.06562  0.16417  1.00561 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.29509    0.10174  -2.900  0.00538 ** 
## no_art_5    -0.02844    0.00529  -5.376 1.67e-06 ***
## c.mal        0.07591    0.19195   0.395  0.69406    
## host1       -0.28841    0.09800  -2.943  0.00478 ** 
## c.mal:host1  0.20087    0.26738   0.751  0.45577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3186 on 54 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.4593, Adjusted R-squared:  0.4193 
## F-statistic: 11.47 on 4 and 54 DF,  p-value: 8.242e-07
an1 <- lm(mas ~ no_art_5, data=n1)
summary(an1)
## 
## Call:
## lm(formula = mas ~ no_art_5, data = n1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63760 -0.20930 -0.06227  0.17180  1.12303 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.37902    0.04887  -7.755 6.00e-12 ***
## no_art_5    -0.03222    0.00368  -8.755 3.75e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3261 on 105 degrees of freedom
##   (39 observations deleted due to missingness)
## Multiple R-squared:  0.422,  Adjusted R-squared:  0.4165 
## F-statistic: 76.66 on 1 and 105 DF,  p-value: 3.748e-14
plot(n1$mas[n1$no_art_5 > 5] ~ n1$no_art_5[n1$no_art_5 > 5])

nfds.stats <- as.data.frame(matrix(nrow = 7, ncol=12))

# Timema abundance
for(i in 1:6) {
  
an.tim <- glm(n1$no_tim ~ n1[[189 + i]] + n1$ln.vol_cub + n1$con_tim,
             family=quasipoisson, data=n1,
             subset=no_tim > 4)
nfds.stats[1, c(2 * i - 1, 2 * i)] <- summary(an.tim)$coefficients[2, c(1, 4)]

}

# arthropod abundance

for(i in 1:6) {
  
an.art <- glm(no_art_5 ~ n1[[189+i]] + n1$ln.vol_cub + n1$c.con_art5*n1$host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
nfds.stats[2, c(2 * i - 1, 2 * i)] <- summary(an.art)$coefficients[2, c(1, 4)]

}

# arthropod richness
contrasts(n1$host) <- c(1, 0)

for(i in 1:6) {

an.rich <- glm(n1$rich_art_5 ~ n1$host * n1[[189+i]] + n1$ln.vol_cub + n1$cn2 + n1$host * n1$c.con_art5, 
           family=quasipoisson, data=n1)
nfds.stats[3:4, c(2 * i - 1, 2 * i)] <- summary(an.rich)$coefficients[c(3, 7), c(1, 4)]

}

# Mass-abundance

contrasts(n1$host) <- c(0, 1)

for(i in 1:6) {
  
an.mas <- lm(n1$mas ~ n1$host * n1[[189 + i]] + n1$host * n1$c.no_art_5 + n1[[189 + i]] * n1$c.no_tim + n1$c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
nfds.stats[5:6, c(2 * i - 1, 2 * i)] <- summary(an.mas)$coefficients[c(3, 7), c(1, 4)]

}

# carbon - nitrogen

for(i in 1:6) {
  
an.cn <- lm(n1$cn2 ~ n1[[189 + i]] + n1$host + n1$ln.vol_cub + n1$no_art_5, data=n1, subset=no_tim > 0)
nfds.stats[7, c(2 * i - 1, 2 * i)] <- summary(an.cn)$coefficients[2, c(1, 4)]
  
}

nfds.stats$term <- c("mal", "mal", "mal", "mal * host", "mal", "mal * host", "mal")
colnames(nfds.stats) <- rep(c("b", "p"), 6)
write.csv(nfds.stats, file = "data/nfds_stats.csv")

Figures

jpeg(file="figures/n1_tim_mal.jpg", 
     width=750, height=750)

par(mar=c(5, 6, 1, 1))

plot(y=jitter(n1$no_tim[n1$no_tim > 4], 2), 
     x=jitter(n1$mal[n1$no_tim > 4], 10), 
     xlab="maladaptation", ylab="", 
     ylim=c(3, 17),
     las=1, pch=16, 
     cex=1.5, cex.lab=1.5, cex.axis=1.5)

title(ylab=expression(paste(italic("T. cristinae"), " abundance")), cex.lab=1.5, line=3.5)

axis.break(axis=2, breakpos=3.5)

dev.off()
## quartz_off_screen 
##                 2
an2.f <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
           family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.f)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host, 
##     family = quasipoisson, data = n1, subset = no_tim > 0)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.6751  -1.6447  -0.7855   0.8684   7.8261  
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.594e+00  9.414e-01  -2.755  0.00714 ** 
## c.mal             4.999e-01  2.369e-01   2.110  0.03770 *  
## ln.vol_cub        4.316e-01  8.042e-02   5.367 6.55e-07 ***
## c.con_art5       -8.682e-05  2.590e-02  -0.003  0.99733    
## host1            -4.782e-01  2.173e-01  -2.200  0.03042 *  
## c.con_art5:host1  6.740e-02  3.701e-02   1.821  0.07201 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
## 
##     Null deviance: 652.94  on 92  degrees of freedom
## Residual deviance: 407.37  on 87  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
an2.f1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
           family=quasipoisson, data=n1, subset = no_tim > 0 & host== "C")
summary(an2.f1)
## 
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson, 
##     data = n1, subset = no_tim > 0 & host == "C")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.6270  -1.5586  -0.7073   0.7373   7.9524  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.05242    1.26538  -2.412   0.0189 *  
## c.mal        0.71362    0.34010   2.098   0.0401 *  
## ln.vol_cub   0.42820    0.09911   4.320 5.95e-05 ***
## c.con_art5   0.07210    0.02852   2.528   0.0141 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 6.078791)
## 
##     Null deviance: 512.65  on 63  degrees of freedom
## Residual deviance: 311.59  on 60  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
n1.2 <- n1[complete.cases(n1$no_art_5, n1$host, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5, n1$no_tim), ]
res.art <- resid(glm(no_art_5 ~ ln.vol_cub + c.con_art5*host,
           family=quasipoisson, data=n1.2, subset = no_tim > 0 ), 'deviance')
res.mal <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5*host,
           data=n1.2, subset = no_tim > 0), "pearson")

an.res <- lm(res.art ~ res.mal); summary(an.res)
## 
## Call:
## lm(formula = res.art ~ res.mal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4703 -1.3778 -0.3566  1.1054  7.9407 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -0.2809     0.2192  -1.281   0.2033  
## res.mal       1.2042     0.6201   1.942   0.0552 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.114 on 91 degrees of freedom
## Multiple R-squared:  0.0398, Adjusted R-squared:  0.02925 
## F-statistic: 3.772 on 1 and 91 DF,  p-value: 0.05522
nd=data.frame(res.mal=seq(min(res.mal), max(res.mal), len=50))
preds <- predict(an.res, newdata=nd, type = "response", se.fit=TRUE)
nd$fit <- preds$fit
nd$lower <- preds$fit - preds$se.fit
nd$upper <- preds$fit + preds$se.fit

jpeg(file="figures/n1_art_mal.jpg", 
     width=750, height=750)

par(mar=c(6, 6, 1, 1))

plot(y=res.art, x=res.mal, 
     xlab="", ylab="", 
     las=1, pch=16, 
     cex=2, cex.lab=2, cex.axis=2)

title(ylab="arthropod abundance (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4.5)
#axis.break(axis=2, breakpos=3.5)

polygon(x = c(nd$res.mal, rev(nd$res.mal)), 
        y = c(nd$upper, rev(nd$lower)), 
        col="grey", density=100, lwd=.5)

lines(nd$res.mal, nd$fit, lwd=2)

dev.off()
## quartz_off_screen 
##                 2
#lines(nd$res.mal, nd$upper, lty=2)
#lines(nd$res.mal, nd$lower, lty=2)
an3.f <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + host*c.con_art5, 
           family=quasipoisson, data=n1)
summary(an3.f)
## 
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + host * 
##     c.con_art5, family = quasipoisson, data = n1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3865  -1.0153  -0.2121   0.7193   2.0200  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.77989    0.63926  -4.349 3.75e-05 ***
## host1             0.02279    0.17607   0.129  0.89730    
## c.mal            -0.25980    0.33514  -0.775  0.44034    
## ln.vol_cub        0.35623    0.05413   6.581 3.51e-09 ***
## c.con_art5       -0.02600    0.02146  -1.212  0.22897    
## host1:c.mal       0.98839    0.39786   2.484  0.01492 *  
## host1:c.con_art5  0.08287    0.02862   2.896  0.00479 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
## 
##     Null deviance: 234.49  on 92  degrees of freedom
## Residual deviance: 122.51  on 86  degrees of freedom
##   (53 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
n1.2 <- n1[complete.cases(n1$rich_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]

# Adenostoma
an3.fA <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, 
           family=quasipoisson, data=n1, subset=host=="A")
summary(an3.fA)
## 
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson, 
##     data = n1, subset = host == "A")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6979  -0.7854  -0.1862   0.5058   1.7825  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -2.76774    1.34874  -2.052  0.05077 . 
## c.mal       -0.25958    0.28588  -0.908  0.37255   
## ln.vol_cub   0.35518    0.11691   3.038  0.00551 **
## c.con_art5  -0.02596    0.01872  -1.387  0.17764   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.9827792)
## 
##     Null deviance: 35.823  on 28  degrees of freedom
## Residual deviance: 24.749  on 25  degrees of freedom
##   (44 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
res.artA <- resid(glm(rich_art_5 ~ ln.vol_cub + c.con_art5,
           family=quasipoisson, data=n1.2, subset = no_tim > 0 & host=="A" ), 'deviance')
res.malA <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5,
           data=n1.2, subset = no_tim > 0 & host=="A"), "working")

an.resA <- lm(res.artA ~ res.malA); summary(an.resA)
## 
## Call:
## lm(formula = res.artA ~ res.malA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6110 -0.7195 -0.1231  0.5514  1.8781 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09272    0.17442  -0.532    0.599
## res.malA    -0.60915    0.46384  -1.313    0.200
## 
## Residual standard error: 0.9393 on 27 degrees of freedom
## Multiple R-squared:  0.06004,    Adjusted R-squared:  0.02523 
## F-statistic: 1.725 on 1 and 27 DF,  p-value: 0.2001
ndA=data.frame(res.malA=seq(min(res.malA), max(res.malA), len=50))
predsA <- predict(an.resA, newdata=ndA, type = "response", se.fit=TRUE)
ndA$fit <- predsA$fit
ndA$lower <- predsA$fit - predsA$se.fit
ndA$upper <- predsA$fit + predsA$se.fit

# Ceanothus
an3.fC <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, 
           family=quasipoisson, data=n1, subset=host=="C")
summary(an3.fC)
## 
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson, 
##     data = n1, subset = host == "C")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3858  -1.1404  -0.2420   0.8255   2.0197  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.75955    0.79599  -3.467  0.00098 ***
## c.mal        0.72862    0.22775   3.199  0.00220 ** 
## ln.vol_cub   0.35643    0.06266   5.688 4.05e-07 ***
## c.con_art5   0.05688    0.01992   2.855  0.00590 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.538553)
## 
##     Null deviance: 184.732  on 63  degrees of freedom
## Residual deviance:  97.756  on 60  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
res.artC <- resid(glm(rich_art_5 ~ ln.vol_cub + c.con_art5,
           family=quasipoisson, data=n1.2, subset = no_tim > 0 & host=="C" ), 'deviance')
res.malC <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5,
           data=n1.2, subset = no_tim > 0 & host=="C"), "working")

an.resC <- lm(res.artC ~ res.malC); summary(an.resC)
## 
## Call:
## lm(formula = res.artC ~ res.malC)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.53433 -0.94271 -0.05687  1.10582  2.21744 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.1589     0.1580  -1.006  0.31827   
## res.malC      1.3378     0.4632   2.888  0.00533 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.264 on 62 degrees of freedom
## Multiple R-squared:  0.1186, Adjusted R-squared:  0.1044 
## F-statistic: 8.342 on 1 and 62 DF,  p-value: 0.005327
ndC=data.frame(res.malC=seq(min(res.malC), max(res.malC), len=50))
predsC <- predict(an.resC, newdata=ndC, type = "response", se.fit=TRUE)
ndC$fit <- predsC$fit
ndC$lower <- predsC$fit - predsC$se.fit
ndC$upper <- predsC$fit + predsC$se.fit

jpeg(file="figures/n1_rich_mal.jpg", 
     width=750, height=750)

par(mar=c(6, 6, 1, 1))

plot(y=c(res.artA, res.artC), x=c(res.malA, res.malC), 
     xlab="", ylab="", 
     ylim=c(min(c(res.artA, res.artC)), 3),
     las=1, type="n", 
     cex=2, cex.lab=2, cex.axis=2)

title(ylab="arthropod richness (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)
#axis.break(axis=2, breakpos=3.5)

points(y=res.artA, x=res.malA, 
       pch=16, col="blue", cex=2)
points(y=res.artC, x=res.malC, 
       pch=16, col="orange", cex=2)

polygon(x = c(ndC$res.malC, rev(ndC$res.malC)), 
        y = c(ndC$upper, rev(ndC$lower)), 
        col="orange", density=100, lwd=.5)

lines(ndC$res.malC, ndC$fit, lwd=2, col="orange")

legend(x=-.5, y=3, 
       legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))), 
       pch=c(16, 16), col=c("orange", "blue"), 
       bg=NULL, bty="n", cex=c(2, 2), y.intersp = 1)

dev.off()
## quartz_off_screen 
##                 2
n1.3 <- n1[complete.cases(n1$mas, n1$c.no_art_5, n1$c.mal), ]

# Adenostoma
an.s2A <- lm(mas ~ c.mal + c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
summary(an.s2A) # R-sq = 18.4 %
## 
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1.3[n1.3$no_art_5 > 
##     2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62088 -0.20001  0.03397  0.15980  0.85752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.40958    0.10117  -4.048 0.000536 ***
## c.mal        0.18082    0.21821   0.829 0.416207    
## c.no_art_5  -0.04825    0.01169  -4.127 0.000443 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3531 on 22 degrees of freedom
## Multiple R-squared:  0.4409, Adjusted R-squared:  0.3901 
## F-statistic: 8.675 on 2 and 22 DF,  p-value: 0.001668
an.mas.artA <- lm(mas ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
an.mal.artA <- lm(c.mal ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])

res.masA <- resid(an.mas.artA)
res.malA <- resid(an.mal.artA)

an.p_malA <- lm(res.masA ~ res.malA)
summary(an.p_malA) # r-sq = 13.5! on Adenostoma
## 
## Call:
## lm(formula = res.masA ~ res.malA)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62088 -0.20001  0.03397  0.15980  0.85752 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.775e-18  6.906e-02   0.000    1.000
## res.malA    1.808e-01  2.134e-01   0.847    0.406
## 
## Residual standard error: 0.3453 on 23 degrees of freedom
## Multiple R-squared:  0.03027,    Adjusted R-squared:  -0.0119 
## F-statistic: 0.7179 on 1 and 23 DF,  p-value: 0.4056
ndA=data.frame(res.malA=seq(min(res.malA), max(res.malA), len=50))
predsA <- predict(an.p_malA, newdata=ndA, type = "response", se.fit=TRUE)
ndA$fit <- predsA$fit
ndA$lower <- predsA$fit - predsA$se.fit
ndA$upper <- predsA$fit + predsA$se.fit

# Ceanothus

an.s2C <- lm(mas ~ c.mal + c.no_art_5 , data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
summary(an.s2C) # R-sq = 18.4 % (8.4)
## 
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1.3[n1.3$no_art_5 > 
##     2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52193 -0.21290 -0.07037  0.18720  0.83286 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.720632   0.047085 -15.305  < 2e-16 ***
## c.mal        0.113283   0.129225   0.877    0.385    
## c.no_art_5  -0.029155   0.004516  -6.456 5.03e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2882 on 48 degrees of freedom
## Multiple R-squared:  0.465,  Adjusted R-squared:  0.4427 
## F-statistic: 20.86 on 2 and 48 DF,  p-value: 3.025e-07
an.mas.artC <- lm(mas ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
an.mal.artC <- lm(c.mal ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])

res.masC <- resid(an.mas.artC)
res.malC <- resid(an.mal.artC)

an.resC <- lm(res.masC ~ res.malC); summary(an.resC)
## 
## Call:
## lm(formula = res.masC ~ res.malC)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52193 -0.21290 -0.07037  0.18720  0.83286 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.999e-17  3.994e-02   0.000     1.00
## res.malC    1.133e-01  1.279e-01   0.886     0.38
## 
## Residual standard error: 0.2852 on 49 degrees of freedom
## Multiple R-squared:  0.01576,    Adjusted R-squared:  -0.004329 
## F-statistic: 0.7845 on 1 and 49 DF,  p-value: 0.3801
ndC=data.frame(res.malC=seq(min(res.malC), max(res.malC), len=50))
predsC <- predict(an.resC, newdata=ndC, type = "response", se.fit=TRUE)
ndC$fit <- predsC$fit
ndC$lower <- predsC$fit - predsC$se.fit
ndC$upper <- predsC$fit + predsC$se.fit

jpeg(file="figures/n1_mas_mal.jpg", 
     width=750, height=750)

par(mar=c(6, 7, 1, 1))

plot(y=c(res.masA, res.masC),
     x=c(res.malA, res.malC),
     xlab = "", 
     ylab = "", 
     las=1, type="n", cex.lab=2, cex.axis=2, cex=2)

title(ylab="mass-abundance slope (residuals)", cex.lab=3, line=4.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)

points(y=res.masC, x=res.malC, 
       pch=16, col="orange", cex=2)
points(y=res.masA, x=res.malA, 
       pch=16, col="blue", cex=2)


#lines(nd$res.mal, nd$upper, lty=2)
#lines(nd$res.mal, nd$lower, lty=2)

polygon(x = c(ndA$res.malA, rev(ndA$res.malA)), 
        y = c(ndA$upper, rev(ndA$lower)), 
        col="blue", density=50, lwd=.5, border=NA)

lines(ndA$res.malA, ndA$fit,
      col="blue", lwd=2)

legend(x=-.5, y=.35, 
       legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))), 
       pch=c(16, 16), col=c("orange", "blue"), 
       bg=NULL, bty="n", cex=2, y.intersp = 1)

dev.off()
## quartz_off_screen 
##                 2
p.n1 <- n1[complete.cases(n1$host, n1$c.mal, n1$cn2) & n1$pN < 1, ]
p.n1$c.host <- ifelse(p.n1$host == "A", -.5, .5)

an.mal <- lm(c.mal ~ c.host, data=p.n1)
an.cn_h <- lm(cn2 ~ c.host, data=p.n1)

r.mal <- resid(an.mal)
r.cn_h <- resid(an.cn_h)

an.res <- lm(r.cn_h ~ r.mal); summary(an.res)
## 
## Call:
## lm(formula = r.cn_h ~ r.mal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3758 -1.3005 -0.0045  1.1843  4.6236 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -3.601e-17  1.950e-01   0.000   1.0000  
## r.mal       -1.301e+00  5.417e-01  -2.403   0.0185 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.798 on 83 degrees of freedom
## Multiple R-squared:  0.06502,    Adjusted R-squared:  0.05376 
## F-statistic: 5.772 on 1 and 83 DF,  p-value: 0.01851
nd=data.frame(r.mal=seq(min(r.mal), max(r.mal), len=50))
preds <- predict(an.res, newdata=nd, type = "response", se.fit=TRUE)
nd$fit <- preds$fit
nd$lower <- preds$fit - preds$se.fit
nd$upper <- preds$fit + preds$se.fit


jpeg(file="figures/n1_cn_mal.jpg", 
     width=750, height=750)

par(mar=c(6, 6, 1, 1))

plot(y=r.cn_h, x=r.mal, 
     xlab="", ylab="", 
     las=1, pch=16, 
     cex=2, cex.lab=3, cex.axis=2)

title(ylab="carbon:nitrogen ratio (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)
#axis.break(axis=2, breakpos=3.5)

polygon(x = c(nd$r.mal, rev(nd$r.mal)), 
        y = c(nd$upper, rev(nd$lower)), 
        col="grey", density=100, lwd=.5)

lines(nd$r.mal, nd$fit, lwd=2)

dev.off()
## quartz_off_screen 
##                 2
# carbon on CN
jpeg(file="figures/n1_cn_pC.jpg", 
     width=750, height=750)

par(mar=c(6, 7, 3, 1))

plot(y=p.n1$cn2, x=p.n1$pC, 
     ylab="", xlab="", 
     las=1, pch=16, col=ifelse(p.n1$host == "A", "blue", "orange"),
     cex=2, cex.lab=3, cex.axis=2)

title(ylab="carbon:nitrogen ratio", line=4.5, cex.lab=3)
title(xlab="mass carbon (mg)", line=4.5, cex.lab=3)

legend(x=7.25, y=33, 
       legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))), 
       pch=c(16, 16), col=c("orange", "blue"), 
       bg=NULL, bty="n", cex=2, y.intersp = 1)

dev.off()
## quartz_off_screen 
##                 2
# nitrogen on CN

an_n <- lm(cn2 ~ pN + I(pN ^ 2), data=p.n1)
summary(an_n)
## 
## Call:
## lm(formula = cn2 ~ pN + I(pN^2), data = p.n1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1372 -0.9234  0.1656  1.0417  4.3042 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   62.181      5.368  11.584  < 2e-16 ***
## pN          -133.353     24.095  -5.535 3.65e-07 ***
## I(pN^2)      100.705     26.418   3.812 0.000266 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.828 on 82 degrees of freedom
## Multiple R-squared:  0.8013, Adjusted R-squared:  0.7964 
## F-statistic: 165.3 on 2 and 82 DF,  p-value: < 2.2e-16
nd <- data.frame(pN=seq(min(p.n1$pN), max(p.n1$pN), len=50))
preds <- predict(an_n, newdata=nd, type="response", se.fit=TRUE)
nd <- cbind(nd, fit=preds$fit, upper=preds$fit + preds$se.fit, lower=preds$fit-preds$se.fit)

jpeg(file="figures/n1_cn_pN.jpg", 
     width=750, height=750)

par(mar=c(6, 7, 1, 1))

plot(y=p.n1$cn2, x=p.n1$pN, 
     ylab="", xlab="", 
     las=1, pch=16, col=ifelse(p.n1$host == "A", "blue", "orange"),
     cex=2, cex.lab=2, cex.axis=2)

title(ylab="carbon:nitrogen ratio", line=4.5, cex.lab=3)
title(xlab="mass nitrogen (mg)", line=4.5, cex.lab=3)

polygon(x = c(nd$pN, rev(nd$pN)), 
        y = c(nd$upper, rev(nd$lower)), 
        col="grey", density=100, lwd=.5)

lines(nd$pN, nd$fit, lwd=2, col="grey")

legend(x=.5, y=30, 
       legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))), 
       pch=c(16, 16), col=c("orange", "blue"), 
       bg=NULL, bty="n", cex=2, y.intersp = 1)

dev.off()
## quartz_off_screen 
##                 2
library(scales)
## Warning: package 'scales' was built under R version 3.6.2
## 
## Attaching package: 'scales'
## The following object is masked from 'package:plotrix':
## 
##     rescale
n1$occ <- ifelse(n1$no_tim_gs==0, FALSE, TRUE)
n1$radius <- (n1$vol_cub/(100000*pi))^(1/3)

n1_occ <- n1[n1$occ,]

adOnly <- n1_occ[n1_occ$host=="A",]
ceOnly <- n1_occ[n1_occ$host=="C",]
empty <- n1[!n1$occ,]

cFunA <- colorRamp(c("skyblue1", "blue"))
cFunC <- colorRamp(c("darkorange2", "peachpuff"))

colsA <- rgb(cFunA(adOnly$mal), maxColorValue=255)
colsC <- rgb(cFunC(ceOnly$mal), maxColorValue=255)

jpeg(file="figures/n1_map.jpg", 
     width=750, height=750)

par(mar=c(1, 1, 1, 1))

plot(x=1:10, y=1:10, 
     xlim=c(0,80), ylim=c(-1, 50), 
     type="n", 
     xlab="", ylab="",
     las=1, cex.axis=2, cex.lab=2, axes=FALSE)
#title(ylab="meters N-S", line=4, cex.lab=2)

polygon(x=c(-3, -3, 85, 85), y=c(-3, 55, 55, -3), density=500, col="lightgrey")
box(which="plot", lty="solid")
polygon(x=c(60, 60, 80, 80), y=c(0, .25, .25, 0), density=500, col="black")
text(70, -.75, labels="10 meters", cex=1.5)
symbols(x=adOnly$x, y=adOnly$y, circles=(adOnly$radius), bg=colsA, inches=FALSE,
        xlim=c(0,75), ylim=c(-1, 50),
        xlab = "metres E-W", ylab="metres N-S", fg="blue", lwd=2.5, cex.lab=1.5,
        cex.axis=1.5, add=TRUE)
symbols(x=ceOnly$x, y=ceOnly$y, circles=(ceOnly$radius), bg=colsC, inches=F, add=T,
        fg="orange", lwd=2.5)
symbols(x=empty$x, y=empty$y, circles=(empty$radius), bg="white",
        fg=alpha(ifelse(empty$host=="A", "blue", "orange"),0.5),
        inches=F, add=T, lwd=2.5)



legend(x=0, y=50, legend=c(expression(italic("A. fasciculatum")), expression(italic("C. spinosus"))), pch=16, col=c("blue", "orange"), bty="n", cex=2)

dev.off()
## quartz_off_screen 
##                 2
ad$taxon <- factor(ifelse(ad$family == "Formicidae", "Formicidae", 
                          ifelse(ad$order == "Homoptera", "Hemiptera",
                                 ifelse(ad$order == "Neuroptera", "Chrysopidae",
                                        ifelse(ad$order == "Aranea",
                                               as.character(ad$family),
                                               ifelse(ad$order == "Coleoptera" & ad$family != "", 
                                                      as.character(ad$family), 
                                               as.character(ad$order)))))))

#ad$taxon <- factor(ifelse(ad$family != "", as.character(ad$family),
#                          as.character(ad$order)))

ad$taxon <- factor(ifelse(ad$family == "Formicidae", "Formicidae", 
                          ifelse(ad$order == "Homoptera", "Hemiptera",
                                 ifelse(ad$order == "Aranea", "Aranea", 
                                 ifelse(ad$order %in% c("Chrysopidae", "Diptera",
                                                        "Hymenoptera"), 
                                        "Other", as.character(ad$order))))))

orders <- levels(droplevels(ad$taxon[ad$length >=5]))

arts <- setNames(lapply(orders, 
                        FUN=function(x) table(ad$plant_id[ad$taxon == x & ad$length >=5])),
                 orders)

newnames <- paste("no_", orders, "_5", sep='')

n1[, newnames] <- 0

for(i in 1:length(newnames)) {
  
  n1[match(names(arts[[i]]), n1$plant_id), newnames[i]] <- arts[[i]]
  
}

pnames <- paste("p_", orders, sep="")

n1[, pnames] <- 0

for(i in 1:length(pnames)) {
  
  n1[, pnames[i]] <- ifelse(n1$no_any_5 == 0, NA, n1[, newnames[i]] / n1[, "no_any_5"])
  
}

thresh <- c(seq(0, .9, length.out = 10), 0.9999999)

surv <- do.call(rbind, lapply(pnames, FUN = function(x) {
  
  sapply(thresh, FUN=function(y) {
    
    sum(n1[, x] > y, na.rm=TRUE)
    
  })
  
})
)

dimnames(surv) <- list(orders, c(thresh[1:10], 1))

arttot <- table(droplevels(ad$taxon[ad$length >= 5]))

surv <- cbind(tot=arttot, surv)

surv <- surv[order(surv[, "tot"], decreasing=TRUE), ]


#ad[ad$order == "Lepidoptera" & ad$length >= 5, c("order", "family")]

#write.csv(surv, file="~/Dropbox/Projects/CA Arthropods/N1/Data/art_surv.csv", 
#           row.names=TRUE)

write.csv(surv, file="data/art_surv_gp.csv", 
           row.names=TRUE)

# smaller things

arttot <- table(droplevels(ad$taxon[ad$length <5]))

ad[ad$family=="Formicidae", c("genus")]
## NULL
library(partitions)

p = 0.05
s = 10

naive.pmf <- function(prob = .05, giveup = 10, max = 40, plot.pmf = TRUE) {

  p = prob
  s = giveup
  
pmf <- c((1 - pgeom(s-1, p)), c(lapply(1:(max-s), function(x) {
  
  library(partitions)
  part <- t(parts(x)) 
  ps <- lapply(
          lapply(1:nrow(part), function(y) part[y,]), function(z) z[z > 0])
  ps <- ps[sapply(ps, function(q) all(q < s))]
  pt.ct <- mapply(FUN=function(unq, len) factorial(len)/prod(factorial(unq)), 
       unq=lapply(ps, table), 
       len=lapply(ps, length))
  probs <- mapply(function(ps, pt.ct) prod(dgeom(ps-1, prob=p)) * pt.ct,
      ps = ps, pt.ct = pt.ct)
  sum(probs) * (1 - pgeom(s-1, p))
  
}), recursive=TRUE))

 if(plot.pmf == TRUE) plot(10:max, pmf, type="l")

cum <- sum(pmf) # cumulative density = 0.948273
expt <- sum(pmf * 10:max) # expected value = 12.0177

list(pmf = pmf, cum = cum, expt = expt)
}

dgeom(0, .05)
## [1] 0.05
dgeom(0, .05)^2 + dgeom(1, .05)
## [1] 0.05
dgeom(0, .05)^3 + 2*(dgeom(1, .05) * dgeom(0, .05)) + dgeom(2, .05)
## [1] 0.05
lm(c(.3, .65) ~ c(.2, .8))
## 
## Call:
## lm(formula = c(0.3, 0.65) ~ c(0.2, 0.8))
## 
## Coefficients:
## (Intercept)  c(0.2, 0.8)  
##      0.1833       0.5833
st.w <- function(freq) { 1.117 - 0.583 * freq }
 gr.w <- function(freq) { 0.183 + 0.583 * freq}

 w.bar <- function(freq) st.w(freq) * freq + gr.w(freq) * (1 - freq)

 wbs <- w.bar(seq(0, 1, .01))

cbind(seq(0,1,.01), wbs)
##                   wbs
##   [1,] 0.00 0.1830000
##   [2,] 0.01 0.1980534
##   [3,] 0.02 0.2128736
##   [4,] 0.03 0.2274606
##   [5,] 0.04 0.2418144
##   [6,] 0.05 0.2559350
##   [7,] 0.06 0.2698224
##   [8,] 0.07 0.2834766
##   [9,] 0.08 0.2968976
##  [10,] 0.09 0.3100854
##  [11,] 0.10 0.3230400
##  [12,] 0.11 0.3357614
##  [13,] 0.12 0.3482496
##  [14,] 0.13 0.3605046
##  [15,] 0.14 0.3725264
##  [16,] 0.15 0.3843150
##  [17,] 0.16 0.3958704
##  [18,] 0.17 0.4071926
##  [19,] 0.18 0.4182816
##  [20,] 0.19 0.4291374
##  [21,] 0.20 0.4397600
##  [22,] 0.21 0.4501494
##  [23,] 0.22 0.4603056
##  [24,] 0.23 0.4702286
##  [25,] 0.24 0.4799184
##  [26,] 0.25 0.4893750
##  [27,] 0.26 0.4985984
##  [28,] 0.27 0.5075886
##  [29,] 0.28 0.5163456
##  [30,] 0.29 0.5248694
##  [31,] 0.30 0.5331600
##  [32,] 0.31 0.5412174
##  [33,] 0.32 0.5490416
##  [34,] 0.33 0.5566326
##  [35,] 0.34 0.5639904
##  [36,] 0.35 0.5711150
##  [37,] 0.36 0.5780064
##  [38,] 0.37 0.5846646
##  [39,] 0.38 0.5910896
##  [40,] 0.39 0.5972814
##  [41,] 0.40 0.6032400
##  [42,] 0.41 0.6089654
##  [43,] 0.42 0.6144576
##  [44,] 0.43 0.6197166
##  [45,] 0.44 0.6247424
##  [46,] 0.45 0.6295350
##  [47,] 0.46 0.6340944
##  [48,] 0.47 0.6384206
##  [49,] 0.48 0.6425136
##  [50,] 0.49 0.6463734
##  [51,] 0.50 0.6500000
##  [52,] 0.51 0.6533934
##  [53,] 0.52 0.6565536
##  [54,] 0.53 0.6594806
##  [55,] 0.54 0.6621744
##  [56,] 0.55 0.6646350
##  [57,] 0.56 0.6668624
##  [58,] 0.57 0.6688566
##  [59,] 0.58 0.6706176
##  [60,] 0.59 0.6721454
##  [61,] 0.60 0.6734400
##  [62,] 0.61 0.6745014
##  [63,] 0.62 0.6753296
##  [64,] 0.63 0.6759246
##  [65,] 0.64 0.6762864
##  [66,] 0.65 0.6764150
##  [67,] 0.66 0.6763104
##  [68,] 0.67 0.6759726
##  [69,] 0.68 0.6754016
##  [70,] 0.69 0.6745974
##  [71,] 0.70 0.6735600
##  [72,] 0.71 0.6722894
##  [73,] 0.72 0.6707856
##  [74,] 0.73 0.6690486
##  [75,] 0.74 0.6670784
##  [76,] 0.75 0.6648750
##  [77,] 0.76 0.6624384
##  [78,] 0.77 0.6597686
##  [79,] 0.78 0.6568656
##  [80,] 0.79 0.6537294
##  [81,] 0.80 0.6503600
##  [82,] 0.81 0.6467574
##  [83,] 0.82 0.6429216
##  [84,] 0.83 0.6388526
##  [85,] 0.84 0.6345504
##  [86,] 0.85 0.6300150
##  [87,] 0.86 0.6252464
##  [88,] 0.87 0.6202446
##  [89,] 0.88 0.6150096
##  [90,] 0.89 0.6095414
##  [91,] 0.90 0.6038400
##  [92,] 0.91 0.5979054
##  [93,] 0.92 0.5917376
##  [94,] 0.93 0.5853366
##  [95,] 0.94 0.5787024
##  [96,] 0.95 0.5718350
##  [97,] 0.96 0.5647344
##  [98,] 0.97 0.5574006
##  [99,] 0.98 0.5498336
## [100,] 0.99 0.5420334
## [101,] 1.00 0.5340000
max(wbs)
## [1] 0.676415
plot(seq(0, 1, .01), wbs)

#wss <- ws(seq(0, 1, .001))
#max(wss - min(wss))
#min(wss)
ws <- function(freq) {
  raw <- 1.517 * freq - 1.166 * freq^2 + 0.183
  (raw - 0.183)/.493415
  }

ws <- function(freq) {
  raw <- w.bar(freq)
  (raw - 0.183)/.493415
  }

ws(seq(0, 1, .01))
##   [1] 0.00000000 0.03050860 0.06054457 0.09010792 0.11919865 0.14781675
##   [7] 0.17596222 0.20363507 0.23083530 0.25756290 0.28381788 0.30960024
##  [13] 0.33490996 0.35974707 0.38411155 0.40800340 0.43142264 0.45436924
##  [19] 0.47684323 0.49884458 0.52037332 0.54142943 0.56201291 0.58212377
##  [25] 0.60176201 0.62092762 0.63962060 0.65784097 0.67558870 0.69286382
##  [31] 0.70966631 0.72599617 0.74185341 0.75723802 0.77215002 0.78658938
##  [37] 0.80055612 0.81405024 0.82707173 0.83962060 0.85169685 0.86330047
##  [43] 0.87443146 0.88508983 0.89527558 0.90498870 0.91422920 0.92299707
##  [49] 0.93129232 0.93911494 0.94646494 0.95334232 0.95974707 0.96567919
##  [55] 0.97113870 0.97612557 0.98063983 0.98468145 0.98825046 0.99134684
##  [61] 0.99397059 0.99612172 0.99780023 0.99900611 0.99973937 1.00000000
##  [67] 0.99978801 0.99910339 0.99794615 0.99631629 0.99421380 0.99163868
##  [73] 0.98859094 0.98507058 0.98107759 0.97661198 0.97167374 0.96626288
##  [79] 0.96037940 0.95402329 0.94719455 0.93989319 0.93211921 0.92387260
##  [85] 0.91515337 0.90596151 0.89629703 0.88615993 0.87555020 0.86446784
##  [91] 0.85291286 0.84088526 0.82838503 0.81541218 0.80196670 0.78804860
##  [97] 0.77365787 0.75879452 0.74345855 0.72764995 0.71136873
plot(seq(0, 1, .01), ws(seq(0, 1, .01)))

freq <- seq(0, .99, .01)
k = 1
ws <- 1- (1 / (1 - exp(freq)))

plot(freq, ws)

freq <- c(0.2, 0.8)
df <- data.frame(a=c(.99, .65), b=I(1 - exp(-freq)))
an1 <- lm(a ~ b, data=df)
preds <- predict(object = an1, newdata = data.frame(b=seq(0, .99, .01)))
k = 1
c = .9

min.RSS <- function(data, par) {
              with(data, sum((par[1] + par[2] * (1 / (1 - c * exp(-x))) - y)^2))
}

dat <- data.frame(x = c(0.2, 0.8), y = c(1, .65))

pars <- optim(par = c(0, 1, 1), min.RSS, data=dat)$par

freqs <- seq(0, .99, .01)
ws <- pars[1] + pars[2] * (1 / (1 - c * exp(-freqs)))
plot(freqs, ws)

nfmal <- function(mal=NULL, k=1) {
  freq <- 1 - mal
  # adapted moph
  min.RSS <- function(data, par) {
              with(data, sum(( (par[1] / (1 - ((par[2] - par[1])/par[2]) * exp((-1) * k * x))) - y)^2))
  }

  dat.ad <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
  pars.ad <- optim(par = c(0,1), min.RSS, data=dat.ad)$par
  
  freqs <- seq(0, 1, .01)
  ws.ad <-  (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freqs)))
  
  # plot(freqs, ws.ad)
  
  # maladapted morph
  
  ws.mal <- rev(ws.ad) - .35 
  
  raw <- (ws.ad * freqs) + (ws.mal * (1 - freqs))
  
  min.w <- min(raw)
  max.w <- max(raw - min.w)
  
  scaled <- (raw - min.w)/max.w
  cbind(freqs, scaled)
  plot(rev(freqs), 1 - scaled)
  #plot(freqs, ws.ad, ylim=c(-2, 2), type="l")
  #plot(freqs, ws.ad, type="l")
  #lines(freqs, ws.mal, lty=2)
  
  # for output
   ws.ad2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freq)))
   ws.mal2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * (1-freq)))) - .35
   raw2 <- (ws.ad2 * freq) + (ws.mal2 * (1 - freq))
   sc2 <- (raw2-min.w) / max.w
  
  #plot(freqs, mal)
  #pars.ad

  1 - sc2

}


#lm(c(1, .65) ~ c(.16, ))
st.w <- function(freq) { 1.117 - 0.583 * freq*(1-freq)}
gr.w <- function(freq) { 0.183 + 0.583 * freq*(1-freq)}

w.bar <- function(freq) st.w(freq-.01) * freq + gr.w(freq-.01) * (1 - (freq-.01))

wbs <- w.bar(seq(0, 1, .01))

cbind(seq(0,1,.01), wbs)
##                   wbs
##   [1,] 0.00 0.1788828
##   [2,] 0.01 0.1941700
##   [3,] 0.02 0.2091085
##   [4,] 0.03 0.2237055
##   [5,] 0.04 0.2379677
##   [6,] 0.05 0.2519024
##   [7,] 0.06 0.2655163
##   [8,] 0.07 0.2788166
##   [9,] 0.08 0.2918103
##  [10,] 0.09 0.3045043
##  [11,] 0.10 0.3169056
##  [12,] 0.11 0.3290213
##  [13,] 0.12 0.3408583
##  [14,] 0.13 0.3524236
##  [15,] 0.14 0.3637242
##  [16,] 0.15 0.3747672
##  [17,] 0.16 0.3855594
##  [18,] 0.17 0.3961080
##  [19,] 0.18 0.4064198
##  [20,] 0.19 0.4165020
##  [21,] 0.20 0.4263615
##  [22,] 0.21 0.4360052
##  [23,] 0.22 0.4454402
##  [24,] 0.23 0.4546735
##  [25,] 0.24 0.4637121
##  [26,] 0.25 0.4725630
##  [27,] 0.26 0.4812331
##  [28,] 0.27 0.4897295
##  [29,] 0.28 0.4980592
##  [30,] 0.29 0.5062291
##  [31,] 0.30 0.5142463
##  [32,] 0.31 0.5221177
##  [33,] 0.32 0.5298504
##  [34,] 0.33 0.5374513
##  [35,] 0.34 0.5449274
##  [36,] 0.35 0.5522858
##  [37,] 0.36 0.5595334
##  [38,] 0.37 0.5666773
##  [39,] 0.38 0.5737243
##  [40,] 0.39 0.5806816
##  [41,] 0.40 0.5875561
##  [42,] 0.41 0.5943548
##  [43,] 0.42 0.6010847
##  [44,] 0.43 0.6077528
##  [45,] 0.44 0.6143661
##  [46,] 0.45 0.6209316
##  [47,] 0.46 0.6274563
##  [48,] 0.47 0.6339472
##  [49,] 0.48 0.6404113
##  [50,] 0.49 0.6468555
##  [51,] 0.50 0.6532869
##  [52,] 0.51 0.6597125
##  [53,] 0.52 0.6661392
##  [54,] 0.53 0.6725742
##  [55,] 0.54 0.6790242
##  [56,] 0.55 0.6854965
##  [57,] 0.56 0.6919978
##  [58,] 0.57 0.6985353
##  [59,] 0.58 0.7051160
##  [60,] 0.59 0.7117468
##  [61,] 0.60 0.7184347
##  [62,] 0.61 0.7251868
##  [63,] 0.62 0.7320100
##  [64,] 0.63 0.7389113
##  [65,] 0.64 0.7458977
##  [66,] 0.65 0.7529763
##  [67,] 0.66 0.7601539
##  [68,] 0.67 0.7674377
##  [69,] 0.68 0.7748345
##  [70,] 0.69 0.7823515
##  [71,] 0.70 0.7899956
##  [72,] 0.71 0.7977737
##  [73,] 0.72 0.8056929
##  [74,] 0.73 0.8137602
##  [75,] 0.74 0.8219826
##  [76,] 0.75 0.8303671
##  [77,] 0.76 0.8389206
##  [78,] 0.77 0.8476502
##  [79,] 0.78 0.8565629
##  [80,] 0.79 0.8656656
##  [81,] 0.80 0.8749654
##  [82,] 0.81 0.8844692
##  [83,] 0.82 0.8941841
##  [84,] 0.83 0.9041170
##  [85,] 0.84 0.9142749
##  [86,] 0.85 0.9246649
##  [87,] 0.86 0.9352939
##  [88,] 0.87 0.9461690
##  [89,] 0.88 0.9572970
##  [90,] 0.89 0.9686851
##  [91,] 0.90 0.9803402
##  [92,] 0.91 0.9922693
##  [93,] 0.92 1.0044794
##  [94,] 0.93 1.0169775
##  [95,] 0.94 1.0297706
##  [96,] 0.95 1.0428657
##  [97,] 0.96 1.0562698
##  [98,] 0.97 1.0699899
##  [99,] 0.98 1.0840330
## [100,] 0.99 1.0984060
## [101,] 1.00 1.1131160
max(wbs)
## [1] 1.113116
plot(seq(0, 1, .01), wbs)

#wss <- ws(seq(0, 1, .001))
#max(wss - min(wss))
#min(wss)
ws <- function(freq) {
  w.bar <- function(freq) st.w(freq-.01) * freq + gr.w(freq-.01) * (1 - freq-.01)
  raw <- w.bar(freq) * freq
  raw/1.10934
  }

ws(seq(0, 1, .01))
##   [1] 0.000000000 0.001717327 0.003701897 0.005944533 0.008436311 0.011168560
##   [7] 0.014132861 0.017321048 0.020725205 0.024337671 0.028151034 0.032158138
##  [13] 0.036352077 0.040726197 0.045274098 0.049989630 0.054866897 0.059900256
##  [19] 0.065084313 0.070413929 0.075884216 0.081490540 0.087228518 0.093094017
##  [25] 0.099083161 0.105192323 0.111418129 0.117757457 0.124207438 0.130765455
##  [31] 0.137429142 0.144196388 0.151065331 0.158034364 0.165102130 0.172267526
##  [37] 0.179529701 0.186888056 0.194342243 0.201892168 0.209537989 0.217280116
##  [43] 0.225119211 0.233056187 0.241092213 0.249228707 0.257467340 0.265810035
##  [49] 0.274258969 0.282816569 0.291485515 0.300268741 0.309169430 0.318191021
##  [55] 0.327337201 0.336611913 0.346019350 0.355563959 0.365250438 0.375083737
##  [61] 0.385069059 0.395211859 0.405517845 0.415992976 0.426643465 0.437475774
##  [67] 0.448496621 0.459712974 0.471132054 0.482761334 0.494608540 0.506681649
##  [73] 0.518988892 0.531538749 0.544339957 0.557401501 0.570732620 0.584342806
##  [79] 0.598241803 0.612439605 0.626946460 0.641772870 0.656929587 0.672427615
##  [85] 0.688278211 0.704492885 0.721083398 0.738061764 0.755440249 0.773231371
##  [91] 0.791447901 0.810102862 0.829209529 0.848781429 0.868832341 0.889376298
##  [97] 0.910427584 0.932000735 0.954110539 0.976772039 1.000000526
plot(seq(0, 1, .01), ws(seq(0, 1, .01)))

freqs <- seq(0, 1, .01)

png(filename = "figures/nfds_mal.png")

plot(freqs, freqs, type="n", 
     xlab= "frequency of cryptic morph", 
     ylab= "average population fitness (scaled)", 
     cex.lab = 1.5, las = 1, cex.axis = 1.5)

for (k in c(0.01, .6 , 1.2, 2, 3, 4)) {
  
  min.RSS <- function(data, par) {
              with(data, sum(( (par[1] / (1 - ((par[2] - par[1])/par[2]) * exp((-1) * k * x))) - y)^2))
  }

  dat.ad <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
  pars.ad <- optim(par = c(0,1), min.RSS, data=dat.ad)$par
  
  freqs <- seq(0, 1, .01)
  ws.ad <-  (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freqs)))
  
  
  # maladapted morph
  
  ws.mal <- rev(ws.ad) - .35 
  
  raw <- (ws.ad * freqs) + (ws.mal * (1 - freqs))
  
  min.w <- min(raw)
  max.w <- max(raw - min.w)
  
  scaled <- (raw - min.w)/max.w
  
  lines(freqs, scaled)
}

text(c(.2, .7), c(.8, .6), 
     labels=c("k = 0.01", "k = 4.0"),
     cex=1.5)

dev.off()
## quartz_off_screen 
##                 2
lm(c(1.667, 1.0835) ~ c(.2, .8))
## 
## Call:
## lm(formula = c(1.667, 1.0835) ~ c(0.2, 0.8))
## 
## Coefficients:
## (Intercept)  c(0.2, 0.8)  
##      1.8615      -0.9725
lm(c(.5, 1.0835) ~ c(.2, .8))
## 
## Call:
## lm(formula = c(0.5, 1.0835) ~ c(0.2, 0.8))
## 
## Coefficients:
## (Intercept)  c(0.2, 0.8)  
##      0.3055       0.9725
st.w <- function(freq) { 1.861 - .972 * (freq/(1-freq)) }
gr.w <- function(freq) { 0.306 + .972 * (freq/(1-freq)) }



w.bar <- function(freq) st.w(freq) * freq + gr.w(freq) * (1 - freq)

wbs <- w.bar(seq(0, .99, .01))
wbs.sc <- (wbs - min(wbs))/max(wbs-min(wbs))

cbind(seq(0,.99,.01), wbs.sc)
##                wbs.sc
##   [1,] 0.00 0.9916345
##   [2,] 0.01 0.9919035
##   [3,] 0.02 0.9921705
##   [4,] 0.03 0.9924352
##   [5,] 0.04 0.9926977
##   [6,] 0.05 0.9929578
##   [7,] 0.06 0.9932155
##   [8,] 0.07 0.9934706
##   [9,] 0.08 0.9937232
##  [10,] 0.09 0.9939732
##  [11,] 0.10 0.9942203
##  [12,] 0.11 0.9944647
##  [13,] 0.12 0.9947060
##  [14,] 0.13 0.9949443
##  [15,] 0.14 0.9951795
##  [16,] 0.15 0.9954114
##  [17,] 0.16 0.9956399
##  [18,] 0.17 0.9958649
##  [19,] 0.18 0.9960863
##  [20,] 0.19 0.9963039
##  [21,] 0.20 0.9965176
##  [22,] 0.21 0.9967272
##  [23,] 0.22 0.9969326
##  [24,] 0.23 0.9971337
##  [25,] 0.24 0.9973302
##  [26,] 0.25 0.9975219
##  [27,] 0.26 0.9977087
##  [28,] 0.27 0.9978904
##  [29,] 0.28 0.9980668
##  [30,] 0.29 0.9982375
##  [31,] 0.30 0.9984025
##  [32,] 0.31 0.9985614
##  [33,] 0.32 0.9987140
##  [34,] 0.33 0.9988600
##  [35,] 0.34 0.9989991
##  [36,] 0.35 0.9991309
##  [37,] 0.36 0.9992552
##  [38,] 0.37 0.9993715
##  [39,] 0.38 0.9994795
##  [40,] 0.39 0.9995788
##  [41,] 0.40 0.9996690
##  [42,] 0.41 0.9997495
##  [43,] 0.42 0.9998199
##  [44,] 0.43 0.9998796
##  [45,] 0.44 0.9999281
##  [46,] 0.45 0.9999648
##  [47,] 0.46 0.9999890
##  [48,] 0.47 1.0000000
##  [49,] 0.48 0.9999970
##  [50,] 0.49 0.9999793
##  [51,] 0.50 0.9999458
##  [52,] 0.51 0.9998958
##  [53,] 0.52 0.9998280
##  [54,] 0.53 0.9997415
##  [55,] 0.54 0.9996349
##  [56,] 0.55 0.9995070
##  [57,] 0.56 0.9993563
##  [58,] 0.57 0.9991811
##  [59,] 0.58 0.9989798
##  [60,] 0.59 0.9987505
##  [61,] 0.60 0.9984909
##  [62,] 0.61 0.9981989
##  [63,] 0.62 0.9978718
##  [64,] 0.63 0.9975069
##  [65,] 0.64 0.9971008
##  [66,] 0.65 0.9966502
##  [67,] 0.66 0.9961511
##  [68,] 0.67 0.9955991
##  [69,] 0.68 0.9949892
##  [70,] 0.69 0.9943158
##  [71,] 0.70 0.9935725
##  [72,] 0.71 0.9927523
##  [73,] 0.72 0.9918467
##  [74,] 0.73 0.9908463
##  [75,] 0.74 0.9897402
##  [76,] 0.75 0.9885157
##  [77,] 0.76 0.9871580
##  [78,] 0.77 0.9856497
##  [79,] 0.78 0.9839702
##  [80,] 0.79 0.9820952
##  [81,] 0.80 0.9799953
##  [82,] 0.81 0.9776350
##  [83,] 0.82 0.9749709
##  [84,] 0.83 0.9719493
##  [85,] 0.84 0.9685033
##  [86,] 0.85 0.9645479
##  [87,] 0.86 0.9599741
##  [88,] 0.87 0.9546390
##  [89,] 0.88 0.9483525
##  [90,] 0.89 0.9408549
##  [91,] 0.90 0.9317830
##  [92,] 0.91 0.9206120
##  [93,] 0.92 0.9065547
##  [94,] 0.93 0.8883742
##  [95,] 0.94 0.8640088
##  [96,] 0.95 0.8297477
##  [97,] 0.96 0.7781690
##  [98,] 0.97 0.6919551
##  [99,] 0.98 0.5191534
## [100,] 0.99 0.0000000
max(wbs)
## [1] 1.088568
plot(seq(0, .99, .01), wbs.sc)

png(filename = "figures/NFDS_model.png")
#par(mfrow=c(2,1))
plot(x=c(0,1), y=c(0,1.2), type="n",
     xlab="cryptic morph frequency",
     ylab="absolute fitness",
     las=1)
abline(a=1.117, b = -0.583)
abline(a=0.183, b = 0.583, lty=2 )
abline(h=.65, lty=3, lwd=.5)
legend("topright", 
       legend=c("cryptic", "conspicuous"),
       lty=c(1,2),
       bty="n", cex = .75)
#plot(seq(0, 1, .01), wbs.sc, type="l", xlab="cryptic morph frequency",ylab="scaled average absolute fitness", las=1)
dev.off()
## quartz_off_screen 
##                 2